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INTRODUCTION 


Work supported under Contract Number NAS8-28248 began in 
January 1972 and has continued with varying degrees of effort to the 
present time. The work has involved the application of improved aero- 
dynamic theory towards the goal of obtaining more accurate knowledge 
of the upper atmosphere. Advancements have been made both in the area 
of aerodynamic theory and the interpretation of the dynamic response 
of objects traveling through the atmosphere. 

The work began as a study of ways of improving models of the 
upper atmosphere as deduced from observation of satellite decay. In 
the development of atmospheric models , it was found that the decay 
of the orbit of a satellite due to drag had been modeled as simply a 
sphere with a drag coefficient of 2.2 traveling through a rotating 
atmosphere. The assumption of a sphere of Cj^ = 2.2 was of course 
recognized to be only an approximation and of particular use in the 
analysis of the drag decay of foreign or some domestic satellites for 
which limited knowledge existed concerning the shape or other physical 
parameters. One of the goals of this work was to investigate the 
magnitude of error made in the assumption of a sphere of C^ = 2.2 and 
to propose more accurate data reduction techniques. 

Chapter 1 describes the major Influence revealed in this study 
concerning the influence of real satellite aerodynamics on the deter- 
mination of upper atmospheric density. Chapter 2 presents a method 
of analysis of satellite drag data which includes the effect of satellite 
lift and the variation in aerodynamic properties around the orbit. 



The method was applied to the data from OVI 15 satellite. One of the 
interesting results obtained from analysis of satellite orbit decay 
has been the super rotation of the atmosphere deduced by King-Hele. 

In Chapter 3, a study is presented which shows that satellite lift 
effects may be responsible for the observed orbit precession rather 
than a super rotation of the upper atmosphere. 

Emphasis of this work gradually came to the lower altitude regime 
and the 80 to 120 Km region in particular. This region of the atmo- 
sphere is of great importance since it serves as a major boundary 
between the lower atmosphere which is constant in molecular weight and 
the upper atmosphere which reaches out many thousands of Km having 
considerable variation in molecular composition. This important region 
of the atmosphere has received little experimental attention due to 
the difficulty and expense of performing measurements at these altitudes. 
The falling sphere method is found to be the primary source of information 
for this region. As with satellite drag measurements, it was found 
that simple assumptions concerning the aerodynamics of objects were 
often employed in the falling sphere analysis. The influence of the 
errors made due to the simplifying assumptions were evaluated and an 
improved method of anlaysis was proposed and applied as reported in 
Chapter 4. 

The work on falling sphere data analysis also revealed that 
most of the 80-120 Km results were based on values of drag coefficient 
in the transition regime that were extrapolated from wind tunnel results 
that were far outside the transition regime. In the work reported here, 
more recent wind tunnel data reported in the literature were obtained 
and more accurate drag coefficient relationships were developed based 
on these data. This work is reported in Chapter 5. 
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The improved drag coefficient relationships revealed a con- 
siderable error in previous falling sphere drag interpretation. These 
data were reanalyzed using the more accurate relationships. The 
results of this work is given in Chapter 6. 

In this work the drag coefficient has been studied for the 
entire spectrum of Knudsen Number and speed ratio. One region which 
was of particular interest is in the very low speed ratio region. 

This region of the aerodynamic spectrum has received little experimental 
attention except for highly viscous flows due to the experimental 
difficulty of obtaining drag data in a low density slow flow situation. 
The theoretical work in this region is discussed in Chapter 7. 

The recommendations for future work are presented in the form 
of two proposals given in Chapter 8. Both proposals would involve 
additional analytic work and subsequent experiments using the shuttle 
space craft. 

The computer programs generated during the period of performance 
of the contract are given in the appendix. 



CHAPTER I 

INFLUENCE OF SATELLITE AERODYNAMICS ON ATMOSPHERIC DENSITY DETERMINATION 

A preprint of a paper by G. R. Karr and R. E. Smith from Preprint 
Volume of the International Conference on Aerospace and Aeronautical 
Meteorology was published by AMS, Boston, Mass. Presentation was given 
May 22-26, 1972. 



■teprinUil from Preprint Volume of the Intimatlonel Conference on Aeroepnee and Aeronsutleel 
Meteorolocr, May 22-20, 1972, Waehlnfton, D.C. ; publMed by AMS, Boeton, Maaa, 


IMFLUfENCE OF SATELLITE AERODYMAMICS ON ATMOS FHEUC DENSITY DETERMINATION’* 


Geral<! R. Karr and Robert E. Smith 

Coordinated Science Laboratory National Aeronavitica and Space Administration 

University o£ Illinois Marshall Space Flight Center 

Urbans, Illinois Huntsville, Alabama 


1 . INTRODUCTION 

Drag'deduced densities of the upper atmos- 
phere have been a primary source of data in the 
development of atmospheric models and to the 
study of the upper atmosphere. In the past, the 
determinstion of atmos^eric density has been 
through the observation of satellite orbital 
decay over a long period of time which neces- 
sarily required knowledge of only the average 
drag properties of the satellite. However, as 
tracking techniques become more accurate and the 
use of sensitive accelerometers increases, the 
assumption of average drag properties la no 
longer valid and a more accurate treatment of 
satellite aerodynamics ssiat be made. 

The purpose of the following discussion will 
be to focus on three principle satellite aero- 
dynamic factors which influence the interpre- 
tation of satellite dynamic response; these are, 
(1) the influence of satellite orientation and 
shape on the drag coefficient, (2) the effect of 
changes in the gas flow properties with altitude, 
and (3) the influence of upper atmospheric winds 
on the interpretation of data. 

The three topics to be treated are effects 
causing the greatest source of error in current 
data reduction. Other factors such as aero- 
dynamic lift, changing atmospheric composition, 
and changing satellite surface properties will 
not be treated here but such factors could be of 
importance for particular satellite systems 
having large aerodynamic lift forces and widely 
varying gas and surface properties. The follow- 
ing will then be limited to a discussion of 
aerodynamic drag effects only and the assumption 
of constant satellite surface properties. The 
atmospheric gas will be considered of single 
species having the average properties associated 
with a particular altitude. 


This work was supported in part by the Joint 
Services Electronics Program (U.S. Army, U.S. 

Navy, and U.S. Air Force) under Contract DAAB- 
07-67-C-0199; and in part by the National 
Aeronautics and Space Administration and the 
American Society for Engineering Education 
through the 1971 ASEE-NASA Suninsr Faculty 
Fellowship Program at Marshall Space Flight 
Center. . ^ 
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The three factors to be discussed are of 
importance in the interpretation of past data as 
well as the application to future more accurate 
measurements. For this reason, an estimate will 
be made of the possible correction to present 
density models based on the results of the 
present study. As a basis for calculation and 
comparison only, the Jacchia (1971) isodel is 
employed in the analysis* Other current models 
could be used for the same purposes with similar 
results. 

2. SATELLITE AERODYMAMICS 

The aerodynamic flow regime experienced by 
Che majority of satellites is free molecular. 

That is, the mean free path between collisions 
of molecules in the upper atmosphere is greater 
than the dimensions of the majority of Earth 
satellites. This assumption is certainly true 
for altitudes greater than 200 km where the near 
free path is on the order of kilometers . 

Although departures from free molecular flow will 
occur in regions of high density, the assumption 
of free molecular flow may be considered 
reasonably accurate to an altitude of 100 Vm 
where the mean free path is of the order of 
meters. For convenience, free molecular flow 
is assumed through out the altitude range being 
considered in this discussion. 

In the free molecular flow of space, by 
definition, collisions of molecules with the 
Satellite surface predominate. For this reason, 
the study of satellite aerodynamics requires an 
understanding of the interaction of gas molecules 
with solid surfaces. The drag properties of a 
satellite are influenced primarily by the 
exchange of momentum with the surface during the 
molecular impact. 

The description of the molecular impact to be 
employed in this discussion is the generalized 
gas surface interaction (CSI) model which has 
direct application to the problem of satellite 
drag, see Karr (1969). This model is represented 
in Figure 1 shewing a general non-specular type 
of reflection. The reflected molecules produce a 
momentum vector in the direction 0. with an 
average velocity of Ui. The reflected properties 
ere assumed to depend'^upon the incident flow 
properties such that 





Uj =■ »/r-aj U 

where and are parameters of Che interaction. 


Generofi^ed 



Fig. I. Description of Che gas surface inter- 
action. 


Using the model of the Interaction described, 
Che force components in Che direction of drag can 
be expressed locally at Che surface. The total 
force is obtained by integrating over the entire 
surface exposed to the flow. If the satellite 
velocity is much higher Chan the random kinetic 
motion of the gas molecules, the gas molecules 
can be considered stationary as the satellite 
sweeps out molecules in its path. The assumption 
of such conditions (Che hypervelocity assumption) 
allows one to determine Che drag coefficient 
based upon the area projected to the flow. 

Under the assumption of hypervelociCy flow 
and Che generalized gas surface interaction Che 
following results are obtained for four shapes of 
interest. The drag coefficient is defined as 

Cjj - Drag/ip V^A 


where A is a reference area* taken to be a 
constant. Independent of angle of attack. 

Flat plate with angle of attack 3 

A “ Area of plate 

Cjj ■ 2 slnB - 2«/1-Q'^ sinP cos j + (2-?^ )3 
Cylinder with axis perpendicular to flow 


diameter 


L ■ length 


[ n n 
cos 2 

(1-Pj)(3-Pj) J 

where the bracketed term is equal to ^ when Pj - 

1.0. A value of P, « 3 is not meaningfully applied 
to this shape . . 



■■ •y 

A “Area of base of cone * n r ■ 


where 6 is the cone half angle. 


Sphere 



where the bracketed term Is equal to zero for 
Pj * 0. A value of “ 4 Is not applicable. 

The above results clearly illustrate Che 
Influence of the GSI and Che shape on Che drag 
properties. The results for a sphere are 
shown plotted in Figure 2 . A range of C- from a 
minimum of 2.0 to a maximum of 4-0 is seen 
depending upon the values of gas surface inter- 
action parameters, and . The results for 

the cone shape reveals that depending upon the 
cone half angle, Cq values less than one arc 
possible. Thus, it Is seen that satellite 
shape and satellite surface properties are 
strong influences on the drag properties. 



Fig. 2. Drag coefficient of sphere as a function 


of gas surface interaction parameters . 

3. THE EFFECT OF SATELLITE ORIENTATION 

Satellites rarely present the same shape Co 
the flow during an orbit. For this reason, the 
effect of satellite orientation, which is 
important in determining the instantaneous drag 
properties of a satellite, is also Important in 
determining the average drag properties. In 
order to illustrate the effect of angle of 
attack, consider first a cylinder with spherical 
end caps with an angle of atcack or. For the 
special case of specular reflection (P. =0) and 
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hyperveloclCy flow, the following analytic 
rcaulta are obtained. 


Cylinder with irh.rtcal ctida ; P. “ 0, hyper* 
velocity flow and angle of attack cr. 


<=0 “ 2 [* ' rl D 

+ -/l -o ^ 7 aln sln^Qr - 2.0*] 

t ^ D Ll J 


_ 2 

where A * ^r and L/D In ihc length to dlaawter 
ratio of Che cylinder. 


These results reveal that C|^,for the cylinder 
with splierlcal ends, changes considerably depend* 
Ing upon the angle of attack and the shape 
(length to diameter ratio). As expected, the Cq 
value Is that of a sphere when a Is zero snd Is 
progressively Influenced by the cylinder as the 
angle of attack is Increased to a value of 

j. Since n is used as the reference area at 

all angles of attack, C|j will reach high valuea 
for large values of L/D- 


For more complex shapes such as the cone, 
results must be obtained numerically. Figure 3 
shows such results for a IS*’ half angle cone. 

The plot Is a three axis presentation of the 
surface Cq(P., 3) where P varies from zero at 
the front to~^a value of *2 at the rear. The 
angle of attack Is varied from zero on the right 
to 180° on the left. Tlie value of a. was taken 
to be zero In constructing this plot. The 
function has a maximum value of 4.00 and a 
mlnlnum value of 0.399. Hiese results Illustrate 
further the strong dependence of aerodynasde 
drag on the shape, orientation, and CSI . 


'o ■ T Ad 

O 

For a circular orbit the satellite will travel at 
a constant rate and, for this case, 

^0 “ I S 

o 

Consider now a spin - stabilized cone shaped 
satellite with a flat base. The orientation of 
the satellite spin axis with respect to the 
orbit Is given by the angle X as shown In 
Figure 4. The angle 6^ serves to orlente the 




Fig. 3. Computer generated plot of surface 

Cjj(Pj,B) for a 35° half angle cone with 

flat base vdiere P. varies from zero to 
two from front to^back and 8 varies 
from zero to IHO’ from right to left. 
The values ol C.. vary from 0.399 to 
4.00. 

Just as the instantaneous values of Cq sre 
Influenced by snglc of attack, the average value 
of Cq over an orbit also depends upon the 
attitude history the satellite experiences over 
the orbit. In order to illustrate this factor, 
consider the determination of the average drag 
coefficient over an orbit of the .satellite. The 
time average of C^ is givui by 
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Fig. 4. Coordinated system describing 

orientation of spin stabilized satellite 
with respect to orbit plane. 


spin axis %rith respect to the velocity vector 
U . For a circular orbit, the velocity vector 
wTll rotate about the spin stabilized satellite 
at a constant rate. Since the instantaneous 
value of Cq for each point in the orbit can be 
determined, the average Cq over the orbit can 
also be determined by numerical quadrature. The 
results of such a study are presented In Figure 
5 for five cone half angles and for X from 0 to 
90°. The base area of the cone was chosen as the 
reference area In each case and Che same CSI 
parameters used for each plot. 

The results given in Figure S Illustrate Che 
Importance of taking the satellite orientation 
history Into account In selecting an average 
drag coefficient. Since the satellite orbit and 
the satellite spin axis will tend to drift In 
space, Che average drag coefficient should not be 
expected to remsln constant. The amount of 
variation Is seen to Increase as Che amount of 
non-symmetry of the satellite Is Increased. Only 
spherically symnetrlc shapes will experience no 
variation. 




Fig. 5. Average drag coefficient for spin 

stabilized cone shaped satellites as a 
function of spin axis orientation with 
respect to orbit. 


A. EFFECT OF ALTITUDE ON DRAG COEFFICIENT (THE 

SPEED RATIO EFFECT) 

The assumption of hypervelocity flow used in 
obtaining the results of the preceding sections 
will now be examined. A discussion of this 
assumption is facilitated by introducing a flow 
parameter termed the speed ratio, S, defined as 
the ratio of the satellite velocity to the 
random velocity of the gas molecules in thermal 
equilibrium. The speed ratio is given by 

S •= U^/y2RT/M 

where R the universal gas constant, M tlie 
av. rage molecular weight and T is the tempera- 
ture of the gas. Hypervelocity flow, the 
assumed flow in the preceding work, is approached 
as the speed ratio approaches Infinity. 

Two factors influence the value of the speed 
ratio as the altitude of the satellite orbit is 
changed. First, considering circular orbits 
only, the velocity of the satellite decreases as 
altitude increases. Second, the temperature of 
the atmosphere and therefore the thermal velocity 
of the molecules increases with increasing 
altitude. These two factors combine to cause a 
considerable decrease in speed ratio as altitude 
increases. This effect is illustrated in 
Figure 6 which shows the value of the circular 
speed ratio as a function of altitude for four 
exospheric temperatures. This plot was 
constructed using values of temperature and 
mean molecular weight from Jacchia (1971) ■ 

In order to investigate the influence of 
speed ratio op the drag coefficient, consider a 
i.-/linder with spherical ends. Utilizing data 
from Karr and Yen (1972) ,Sentman (1961), and 
Fan and Andrews (1969), the value nf for 
throe speed ratios were obtained giving the 
following results for the zero angle of attack 



Fig^ 6. Circular speed ratio as a function of 
altitude for four exospheric tempera- 
tures based on Jacchia 1971 model. 


case. These three values were used to obtain a 



® 2.0 


8 2.106 + 0.450(2l/ttd) 

4 2.249 + 0.885 (2 L/nD) 

second order polynomial approximation to the 
variation In Cp with S. The results are 

Cp/C » 1 + (0.350 + 1.166 L/D)/S 

,+ (0.592 - .1528 L/D)/S . 

This function Is plotted in Figure 7 for three 
values of L/D. Two factors of importance arc 
illustrated in these results. First , long 
slender object at low angles of attack ace 
strongly Influenced by speed ratio effects. 

Even at relatively high values of speed ratio, 
there is strong sensitivity to the length to 
diameter ratio for long slender objects. 

Second , for a given satellite shape, consider- 
able change in the value of Cp is seen to occur 
over the range of interest from S=4 to S = 20. 

Combining the results of the polynomial fit 
of Cp values for a sphere with the change of S 
with respect to altitude, the change in Cp for a 
sphere over the altitude range from 100 to 1000 
km Is obtained and presented in Figure 8, These 
results, for the same four exospheric 
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Fig. 7. Change In drag coefficient as a function 
of speed ratio for three values of L/D. 



Fig. 8. Change in drag coefficient as a function 
of altitude for four exospheric 
temperatures. 

temporatuc'is as in Figure 6, show that a sphere 
will experience a change in Cp of from 2 to 10% 
depending upon the atmospheric temperature. 

Even greater changes would be expected for non- 
spherical satellites such as the cylinder with 
spherical ends . 

Of importance is the fact that the noted 
change in with altitude is found to be 
systematic with altitude. For this reason then, 
it is to be expected that present drag deduced 
density values will have these systematic errors 
incorporated into the results. This factor 
could account for some of the discrepancies in 
densities at high altitude in comparison to 
densities at lower altitudes. Similar results 
and conclusions were found by Izakov( 1965) 
using an analytic expression for Cp of a sphere 
as a function of S. 


S. EFFECT OF ATMOSPHERIC HINDS ON DRAG STUDIES 

Recent analytic and experimental studies 
point to high velocity atsK>spheric winds in the 
upper atmosphere. Velocities of as much as 
AOO m/sec have been reported. Since satellite 
velocities are of the order 7 km/aec, atmospheric 
winds are expected to be an influence on Che 
satellite motion. In order to investigate such 
influence, consider a simple model of the upper 
ati':;Ospheric wind structure which includes the 
rotation of the atmosphere at the rotation rate 
of the Earth. In addition to the rotation, 
consider an east-west wind which varies in both 
magnitude and direction. Results of Challinor 
(1969), suggest that the east-west component 
could be assumed to approximate a sinusoid 
variation with a longitude angle, a, measured 
from a reference point of zero wind. The 
velocity of a Satellite with respect Co the 
atmospheric gas in circular orbit is then given 
by 

U “ ~ rCi - V sinOf 

“ re max 

where V is the peak wind velocity, r is the 
distance from the center of the earth and 0^ 
is the angular velocity of the earth. At some 
longitudes, Che wind is seen to subtract from 
the atmospheric velocity while adding at other 
longitudes. 

The instantaneous drag acting on a satellite 
at any point in the orbit is given by 

D - i PU^C^A. 

The average drag over one revolution is found by 
Integrating over the time in the circular orbit 

* i PC_A I + 5 1 

D L r e max I 

the average drag is found to be increased by a 
factor, where 



r "e 


The increase in drag due to winds for a circular 
orbit is found to be small enough to be 
neglected. The Increase is less than 1% for 
typical values of velocities. 

Consider now a more severe case when a 
satellite in an elliptic orbit has its perigee 
at the peaks of wind velocity. Such a situation 
is likely since, often times, observations of 
elliptic satellite orbits are made to determine 
atmospheric properties in the region of perigee. 
Much of the knowledge of latitude and longitude 
variations in the atmosphere have developed 
from such observations. 

Consider the perigee passage at O' = 90°, 
where the wind velocity subtracts, with the 
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perigee passege et 270°, where the wind velocity 
adds. For almllar orblta and constant denalty, 
Che ratio of drag at perigee for theae two eaaea 
could be taken Co be 



where r la Che orbit radluG at perigee ai\d e la 
the eccentricity of the orbit. In Che altitude 
range from 100 to 1000 km, r^l^^ la approximately 

.5 km/aec, and </n/rp la approximately 7.5 km/aec. 
For e valuea leaa than 0.2, percentage 
dlfferencea In drag of the order of three per- 
cent are obtained for V-^ y valuea of 200 to 
400 m/aec. The wind ef»ct la found to be much 
larger for thla caae than for the circular orbit 
case. Due to Che lag In aolar heating of Che 
atmoaphere, the u-90° point would occur after 
sunset while the 270° point would be after sun- 
rise. The results obtained reveal that a 
difference In drag of a maximum of 2i% at about 
200 km la explainable by wind effects. If 
these wind effects were not taken Into consldera- 
Clonj the difference In drag would be misinter- 
preted aa being caused by a corresponding 
difference In atmospheric density. Therefore, 
wind effecCB could explain some of the day-night 
variation in denalty deduced from aatelllte 
drag. 


6. CONCLUSIONS AND DISCUSSION OF RESULTS 

The preceding discussion has emphasized the 
variability of satellite drag coefficients. In 
particular, the effect of the gas surface 
Interaction Is seen to dominate. Unfortunately, 
values of P. and a have yet to be determined 
accurately enough to be used In satellite drag 
studies. This factor leads to considerable 
uncertainty In the specification of aatelllte 
drag coefficients which could cause errors of 
as much of 507. In the values of Cq “ 2.2 used 
In past data reductions. 

In view of the Influence on drag coefficients 
of non -spherically syimnetrlci shapes and the 
speed ratio, values of of: around 2. '2 are 
likely too low for most satellite shapes. This 
observation is based on the fact that for a 
Sphere the minimum Cq for s-** la 2.0. Speed 
ratio effects cause Che minimum value to be 
Increased to values of 2.1 or higher. In order 
for Che sphere to have the ilnlmuin Cq, Che gas 
surface Interaction would have to be speculat 
or the accommodation coefficient would have to 
be unity corresponding to reflection at near 
zero velocity. Such limiting Interactions 
appear unlikely, meaning that Cq values higher 
than 2.2 are likely. 

Although an absolute value for Cq may be 
lacking presently, the results obtained for 
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changes of Cq due to shape, orientation, and 
altitude reveal that ayetematlc variation of Cq 
around an absolute value would occur under most 
elrcuaatanees. These variations should be 
Included In the analyele of future satellite 
dete where Information on satellite orientation 
may be available. 

The effects of atmospheric winds were 
Illustrated assuming a simple R»del of Che wind 
structure. These results Illustrate still 
mrather source of error In drag deduced 
density values. Future planned work will 
consider a more sophisticated model of the wind 
structure. This study Is expected to lead to 
Improved knowledge of the upper atmospheric 
denalty and wind structure. 
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SATELLITE DYNAMIC RESPONSE 



SATELLITE AERODYNAMICS AND ATMOSPHERIC DENSITY DETERiMlWTI 0 N FROM 

SATELLITE DYNAMIC RESPONSE 

by 

Gora) d R. Karr- 


ABSTRA CT 

A method for determining satei 1 i te acrodynami c properties and upper 
atmospheric density from observed sate I 1 |te dynamic response has been 
successfully developed and tested. 

The aerodynamic drag and lift properties of a satellite are first 
expressed as a function of two parameters associated with gas-surface 
interaction at the satellite surface. The dynamic response of the 
satellite as it passes through the atmosphere is then expressed as 
a function of the tv/o gas-sui'face interaction parameters , the 
atmospheric density, the satellite velocity, and the satellite 
orientation to the high speed flow. By proper correlation of the 
observed dynamic response v»fith the changing angle of attack of the 
satellite, it is found that the tvjo unknown gas-surface interaction 
parameters can be determined. Once the gas-surface interaction 
parameters are known, the aerodynamic properties of the satellite at 
all angles of attack are also determined. The atmospheric density 
may then be accurately calculated once the true aerodynamic properties 
are known. 

Employing accelerometer data from the 0Vl-l|j satel-lite, analysis was 
successfully made of the aerodynamic properties of that satellite 
and a determination was made of the absolute value of atmospheric 
density near the orbit perigee. These results constitute the first 
successful application of the proposed method of analysis. These 
results also serve to illustrate the potential of the technigue in 
the analysis and prediction of satei I i te orbi t ..decay in the atmos- 
phere and the accurate determination of upper atmospheric density 
from satei 1 i to dynam ic response. 




1 . Introduction ' 

The problem of satellite orbit decay prediction and the problem of 
upper atmospheric density determination have encountered a common 
source of unknown error which can be traced to a lack of knowledge of 
satellite aerodynamics. The basic equation employed in both orbit 
decay and density detei-mination is the familiar drag equation: 

Drag = 1/2 pU^CpA 

where p is the density, U is the velocity of the satellite with respect 
to the atmosphere, Cj) is the drag coefficient, and A is a suitable 
reference area. In most applications of this equation to satellites 
the value of Cj^ is considered to have a constant value. Generally, 
however, the assumption of constant drag coefficient is not valid and 
the use of such an assumption ca.n lead to considerable error (see 
Karr, 1972). The uncertainty in satellite aerodynamics has prevented 
the assignment of even an approximate value of drag coefficient with a 
known range of uncertainty, A value of Cp) of 2.0 or 2.2 is often used 
in satellite drag studies and in the determination of atmospheric density. 
These values of Cq are likely too small and, combined with the fact 
that Cp) is not constant, have resulted in an overestimation of upper 
atmospheric densities (see Karr and Smith 1972). 

A more accurate treatment of satellite aerodynamics has obvious 
benefit to the determination of upper atmospheric density and the pre- 
diction of satellite orbit decay. Satellites traveling in the earths upper 
atmosphere experience the aerodynamic flow regime termed the free 
molecular flow regime. In this flow regime, the collision of atmospheric 
gas molecules with the satellite surface dominate the flow and collisions 
of gas molecules with other gas molecules may be ignored. Satellite 
aerodynamic properties are then the result of the interaction of high 
speed gas uiolecules with the solid satellite surface. Unfortunately, 
very little is known about the gas surface interaction at satellite veloci- 
ties and this lack of information is the basic source of uncertainty in 
satellite aerodynamic properties, . . 

In the interest of developing a more accurate treatment of satellite 
aerodynamics for application to orbit decay prediction and density deter in 
nation, a model of the gas surface interaction has been developed which 
utilizes two parameters to describe the interaction (see Karr, 1969 and 
Karr and Yen, 1970). The advantage in this treatment of satellite 
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aerodynamics is that no a priori assumptions of the aerodynamic 
properties need be made. The gas surface interaction parameters 
are considered as unknowns to be determined from the observed 
dynamic response of the satellite as it travels through the atmosphere. 
The determination of the gas surface interaction parameters serves 
as the key to the subsequent determination of both the aerodynamic 
properties and the atmospheric density. 

To tost the proposed method of analysis, accelerometer data from 
the OVI-15 satellite is used. The accelerometer data provide an 
accurate, instantaneous measure of the level of aerodynamic force 
and the attitude of the satellite with respect to the flow. As is pointed 
out in the paper, accurate satellite attitude infor matin is essential to 
the analysis. The OVI-15 satellite, although less than ideal in shape 
for an aerodynamic study, provided a good basis for the test of the 
proposed method of analysis. The results serve to illustrate the 
potential that this method of analysis has to future determinations of 
aerodynamic and atmospheric properties. 

2. Satellite aerodynamics and the gas surface interaction. 

Consider a local satellite surface element in which the high speed 
flow of molecules is incident at an angle of 0 as shown in Figure 1. 
Associated with the incident flow is the incident momentum which 
gives rise to the incident force This force is colinear with the 

satellite velocity, U, with respect to the atmosphere. Assume for 
now that the speed ratio is infinite where the speed ratio is defined 
as the satellite velocity divided by the thermal velocity of the gas 
molecule s. The thermal velocity of the gas molecules is taken to be 
equal to R T/M where R is the gas constant, T is the temperature, 

M is tlie mean molecular weight. 

The molecules reflected from the surface cause a net reaction 

... 

force F j. which is colinear with the mass-motion velocity vector Uj 
of the molecules leaving the surface. The direction of Uj is given 

by the angle G •. 

J 

Modeling of the interaction is performed by providing relationships 
between the incident and reflected quantities. The relationships are 
given by 

U . ■- V 1 -a. U 

^ I’j + (1-Pj) 0 
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where 0!j and Pj are the parameters of the interaction. The subscript j 
is used if more than one set of interaction coefficients is to be considered 
in the analysis. Until more is learned of the interaction, these linear 
relationships provide a useful first approximation to the interaction that 
occurs at satellite velocities. The parameters CZj and Pj are capable 
of describing a much wider range of possible interactions than other 
models. The development of this model and the capabilities are 
described in detail in Karr 1969 and Karr and Yen, 1970. 

The niodcl described above is particularly useful in the determini- 
nation of forces acting on the satellite surface. The total vector force 
acting on the element of surface shown in Figure 1 is given by 

dF = -(U-ajUj) p U-n dA 

where a j is employed if more than one gas surface interaction is 
employed in the analysis. In order to conserve mass at the surface, 
the sum of the o j values must be unity. 

The magnitude and direction of Uj is determined by the parameters 
aj and Pj. The vector force acting on the local element of surface is 
then expressed as a function of 0£j, Pj, U, p, and 0 . For a given satellite 
shape (assumed to be convex) the total aerodynamic forces and torque 
acting on the satellite are found by integration of dl? and RxdF over the 
surface exposed to the flow. In general, the results will be of the form 

Drag = 1/2 p Cp (0!j. Pj,P) A 

Lifti 2 = 1/2P U2 Clj 2<“j' 

Torquej^2,3 = Cxj ^ 

where fi is an angle of orientation and the subscripts on Cj^and are 
to indicate that there are two components of lift and three components 
of torque. The six cierodynarnic properties are found to be a strong 
function of the gas surface interaction parameters. For non-spherical 
objects, the angle of orientation, P , .tlso has a strong influence on the 
drag, lift and torque properties (sec Karr and Yen, 1970). 
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3, Aerodynamics of the OVI-15 satellite. 

The approximated shape of the OVT-15 satellite is a cylinder with 
spherical ends as shown in Figure 2. The satellite spin axis was 
normal to the longitudinal axis of the cylinder with a spin rate of about 
10 rpm. Near the center of the satellite a three axis accelerometer 
detected the forces acting on the satellite. Since the data to be used in 
the subsequent analysis lias been filtered and averaged over a number 
of spin cycles^the aerodynamic properties averaged over a spin cycle 
are developed. 


by 


The total instantaneous vector force acting on the satellite is given 


?• = 1/2 P (Cd’d + Clj Li + Cl2 ^2) 


where Di Lp and are mutually orthogonal unit vectors in the drag 
and lift directions. The Lj and directions are defined v/ith respect 
to the instantaneous orientation such that L] is perpendicular to the 
cylinder axis and the velocity vector. The direction of L 2 is perpendi- 
cular to bath Lj and D. Due to symmetry the lift force in the Ll 
direction is zero. 


From Karr, 1969, Cp and Cl 2 obtained for the infinite speed 
ratio case, given by 

Cj 3 = 2 + 4 Vl -aj (1 -cos Pj (4 - pj) 

+ 2 Ap cos 0s 

+ V 1 -a; J OsSln? Cj - co»3 0. 

R J IT 

sin^ i (Cj + Sj) ]d I 

2it 



where the first two terms in Cj^ are due to the spherical ends and the 
remaining terms are due to the cylindrical section. The angle § is a 
cylindrical surface-integration angle. The quantity Ap is the area ratio 
of the cylinder to the sphere given by 

Ap = ^cyl/-^sph ~ 2rL/*^r = 4 l/tt D 
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The quantities C: and S- contain the parameter P? where 
J J J 


cos (u/2 Pj -f (1 - PQ 9) 




cos 9 

slnC'/2 Pj ^ (• - Pj) e) 


sin 0 


0 = sin”^ (— C03 0g sin?) 


The angle 0 is an angle of instantaneous orientation defined as the angle 
between the velocity vector and the longitudinal axis of the cylinder. 

Since the OVI-J5 spin axis is perpendicular to the axis of the cylinder, 
the angle 0 s function of the angle the spin cucis makes with the velocity 
vector, Y, and a spin angle, X , which changes from 0 to 2rr every spin 
cycle (see Figure 3). 

For certain values of Pj the surface integrals over the angle 5 are 
easily performed. For Pj = 0, which corresponds to specular type 
reflection. 


= S, 


Pj=0 


= 1 


Pj=0 


For P- -- 1 which corresponds to diffusive type reflection, 

f} 




= 0 


s. 


Pj=l 


P:=l 


= -l/cos 0 sin ? 


For P: = 2 which corresponds to perfect backscatter type reflections 


= -I 


Pj=2 


S.| = 1 

j|Pj=2 


The values of and Cj^2 three values of P- = 0, 1, 2, were 

used to obtain an polynomial approximation for Cj) and Ct as a function 
of Pj. 

since the accelerometers were body fixed, the outpxit of the acceler- 
ometers were a function of both drag and lift forces given by 


f/1 pu^ a = 


-Cj^siny cos X + Cj 


1 _ sin^ Y cos ^ X 
cos 0g 
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2 

^ _ 'sin Y cos X sin X 

- Cpj Bln Y sin X - C, 

° L cos 

[ sin Y cos Y ^ 

- C£) cos Y ” CjJ cos 6 ^ 

where i, j, and k are unit vectors in the directions of the three axis of 
the accelerometers with i along the axis of the cylinder, k is in the 
direction of the nominal spin axis and j is orthogonal. Since the data 
being used in this analysis has been averaged over spin cycles, the 
component of force as expressed above were integrated over the angle X. 
The results after averaging oyer one spin cycle were 

Fy= sih'Y^ PU^ nr^c^iOf., Pj, Y) 

= COS Y \ PU^ TTr^ Cp («j, Pj. Y ) 

where is the integrated force coefficient. * These results show that 
Fy and F^ measure the identical forces except for the factor sin y 
cos y, TViis property was used by Fess and Young to obtain the angle Y 
which the spin axis makes with the velocity vector. 

Y = cos-’ /7r/ + f/] 

The force coefficient Cp' v/as found by fitting a 3^^ order polynomial 
to the values of Cp and Cl 2 three values of Pj = 0, 1, and 2. 

2 3 

Cp = A 4v'"l - aj (G+HPj + QPj P Pj ) 

where A = -2 -4 Aj^ E(Y,^) /n 

G= -4Aj^ E (Y. ^)/3tt 

H = 4 F - 2G - 2A 

Q = -4 F 4 5g/4 4 11 a/4 

P = F - Q/4 - 3 a/4 

F“"4/3 - tiAj^/ 2 
n 

where E(Y, -^) is a complete clleptic integral of the second kind resulting 
from the average over one spin cycle. 

^ The accelerometer in the x direction did not function so only fy and 

P are treated in the analysis, 
z ’ 
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4, Method of anaJysis,’ 


The objective of the analysis is to find the best values ofQJj, Pj and 
density which explains the observed accelerometer output of the OVI-15. 
Although Ofj, Pj, and P are considered as unknowns in the analysis, 
certain assumptions on the characteristic variation of p are made to 
facilitate the analysis. The assumption made is that the density varia- 
tion is symmetric with respect to the perigee of the orbit. The absolute 
value of density is still treated as an unlcnown quantity. 

4. 1 Least squares fit. 

Assuming symmetrical density variation about perigee, differ- 
ences in forces measured at points equal distance from perigee must 
be due to changes in the aerodynamic force coefficient, Cp. Since the 
density iff equal at these two points, we can write 

^ P, (T-Atj) U2^02(T +At.) 


where the subscript 1 indicates approach to perigee, subscript 2 indicates 
recession from perigee and T is the perigee passage time. The aero-- 
dynamic properties and forces measured at these two points must then 
satisfy the following relationship 

j^il 

CF(orj. Pj, Yii)’ Cp {ay Pj. Yiz) 

where Yu Yi2 are the angles of orientation at T-«,Atj and T +At^ 

respectively, and h =y 1'.^-^ ^ y^ ' analysis, the quantity DEL^ is 

found from the proceeding relation, defined as. 


DEL; = Til Cp.^ 




il 


where i is used to indicate a comparison made at T +. Atj, A solution in 
the least squares sense is obtained by finding the values of etj and Pj 
which provide a minimum to the sum-of-DELj- squared for a number of 
observations near perigee 

n 

SUM ^ X, (DEL,)2 

1=1 t' 


The best values ofex; and P- are those which satisfy 

J J 

„ ^(SUM) 

in the region of 0 ^ Pj 2, and 0 ^ “f^j ^ 2. 
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4. 2 Perigee passage time. 

The analysis requires first that the aerodynamic properties be 
different at the comparison points used in the analysis. If the Cp were 
not different, then the last equations would be satisfied for all vjuues of 
ai and Pj. Although the OVl-15 satellite was designed to maintain a Y = 

9o° thoughout the orbit, considerable uncontrolled drift in the spin axis 
was found to occur. This malfunction was desirable for purposes of this 
analysis since the angle Y changed considerably during a given orbit. 

This factor resulted in a changing value of Cp over the orbit which pro- 
vided a good sampling of Cp and F values for a wide range of angles of 
orientation. The analysis also requires that the perigee passage time 
be accurately known since the comparisons are made at equal points on 
each side of this time. Since tlie report of Fess and Young did no.‘ provide 
a perigee passage time, it was necessary to calculate that time fr'^m the 
accelerometer output. Since the aerodynamic properties are changing 
during a perigee pass due to the changing angle of orientation, the peak 
in the F curve is shifted in time from the peak in dynamic pressure. 

The derivative of F during a perigee pass is 

- ' II 

F = S Cp + S Cp 

where S is the dynamic pressure. A^t the perigee passage time, the 
dynamic pressure is maximum and S = 0. Therefore, at the perigee 
passage time 

F(T) 

F (T) = Cp-(T) Cj, (T) 

This equation was employed to find the true value of T for each data set 
employed. The quantities Cp and Cp are a function ofOf . and Pj in 
addition to the angle Y. The analysis was able to take info consideration 
the expected shift in F output which was found to vary from 3 to 15 seconds 
depending upon Ofj, Pj and the rate of change in the angle ^ . 

4. 3 Speed ratio effect. 

As discussed by Karr 197 2 and Karr and Smith 197 2, changes 
in speed ratio with altitude result in a systematic increase of Cd with 
altitude. The amount of increase was found to be a function of the 
satellite shape and orientation. This factor was taken into consideration 
in the analysis of data of the OVI-15 by approximating the expected cliange 
in aerodynamics properties with respect to speed ratio. Using information 
from Karr and Smith 197 2 and taking into account the average over a spin 
cycle, the following speed ratio correction factor was obtained. 

COF = 1 4 . 682/S+ . 56148/S^ 

4 . 4 sin2 Y ( 1 . 66/s - . 1528/s^) 
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The correction factor is a function of the angle of orientation which must 
be taken into account in the determination of Ctj and p.. For the altitudes 
of interest for the OVI-15, speed ratios of about 10 are obtained which 
result in an increase in the aerodynamic force coefficient of about 5%. 

The determination of density is then strongly influenced by the speed 
ratio effect. 

5, Results. 

5. 1 Data used in the analysis. 

An example of the data employed in the analysis is shown in 
Figure 4. Data from orbits number 890, 893 and 896 were employed in 
the least squares fit to CKj and Pj. These orbits had significant changes 
in Y during perigee passage and experienced approximately the same 
atmospheric conditions. These nearly polar orbits occurred in mid- 
September 1969 with perigees at 150 km, perigee latitude at 10 ®S latitude 
with a local perigee time at about 1930. Since the sun declination 
was near +3*, the variation in density with latitude near perigee is 
approximately zero according to the Jacchia 1971 model of the atmosphere. 

The choice of orbits used in the analysis contributed to the 
reduction of errors resulting from any nonsymmetry of density variation. 
Further reduction in this type error was made by using only data within 
±10" of perigee. In addition, since data from three orbits was used, 
the errors due to wave motion or other short duration density disturbances 
would contribute only to the random error. 

The values of accelerometer output P, angle of orientation Y , 
and time in seconds from the beginning of data transmission are given 
in Table 1. The units on F is in counts in which 5. 3 counts equals 10”°g 
of acceleration. The angle ^ is given in degrees. The data is seen to 
cover an angle of attack range of about 28 degrees. All the data falls 
within 150 seconds of perigee which for these orbits means that the data 
is taken within ± 10 ® of perigee. Since the true perigee is always within 
± 20 seconds of the peak F, corrected perigee times will not change the 
range of data significantly. 

5. 2 Gas Surface interaction parameters. 

Using the data given in Table 1 and taking into account the perigee 
passage time correction and the speed ratio correction, the results of 
the sum-of-the-squares -of-DELi are given in Figure 5. These results 
show a uni que minimum of the sum-of-the-squares -of-DEL^ at Pj = . 44 
and / 1 -O^j = . 6. These values for Ofj and P; mean that the reflection is 
between a specular and diffusive type reflection in direction and has 
moderate accomadation of energy. 
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5. 3 Aerodynamic properties. 


Using eXj = .64 and Pj - .44, the aerodynamic force properties 
of the OVI-15 satellite may be found using the equations already derived. 
The results are shoNvn in Figure 6 which gives Cp as a function of angle 
of orientation. The level of Cp is dependent directly on the reference 
area, K, choosen to represent the satellite. The plot is given for two 
acceptable areas (1) the maximum cross sectional area seen by the flov^ 
irr^+ZrL and (2) the minimum cross sectional area seen by the flow, 
These results are for the infinite speed ratio case and must be modified 
according to speed ratio influence. 


At angles of 0, 90 and 180°, the quantity Cp is equal to the drag 
coefficient. At all other angles, Cp is influenced by both drag and lift 
coefficients. These results sliow that the drag coefficient is higher than 
the value of 2.2 normally assumed in drag analysis. The speed ratio 
correction will cause these values to be increased by about 5% to 10% for 
the altitudes at which the data was taken. 


The atmospheric density values may be obtained since the values 
of Cp throughout the orbit have been calculated. Assuming an orbit of 
c = .113 and perigee at 150 km, the value of U at the data points arc 
obtained and density values are given by 

P(I, J) = 2am F(I, J)/U^(I, J)A‘Cf(1, J) COF 

where COF is the speed ratio correction factor dependent upon the angle 
Y(I, J), a is the conversion factor needed to convert accelerometer counts 
into accelerometer values (a = lO^g/5.3 counts), and m is the satellite 
mass = 214 kg. The speed ratio correction requires an estimate of the 

speed ratio. For the date, time, and region of the atmosphere for which 
the data corresponds, an estimate was made of the exospheric temperature 
from information provided by Smith, 197 2, In this region of the atmosphere 
the exospheric temperature remains essentially constant and was tciken to 
be 1100®. Using the Jacchia 1971 model atmosphere, a value of t/m versus 
altitude were fitted to a polynomial over the altitude range of interest from 
140 km to 220 km. The speed ratio is then given at each data point by 

S = V/J 2R t/m . 

where K is the universal gas constant. Values of density calculated in 
this manner arc about 5 to 15% less than those predicted by the Jacchia 
1971 model. A more complete discussion of these results will be made 
at a later date. 

6. Discussion of results. 

The results of this analysis are intportant for a number of reasons. 

the analysis illustrates a new method for the analysis of satellite. 
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dynamic response. Second, llic results obtained forOfj and Pj are the 
most accurate of measurements of the gas surface interaction parameters 
at satellite velocity. Third, the values of Cp obtained in the analysis 
are the most accurate of measurements of satellite aerodynamic proper- 
ties. Finally, the results obtained for density are the most accurate 
of measurements of absolute values of upper atmospheric density. 

Past drag analysis have required that critical assumptions be made 
on or p or some of the gas surface interaction parameters. The 
analysis presented here on the other liand has employed very few assump- 
tions in comparison. The assumption of symmetrical density variation 
about perigee is most subject to error. However, the possible error 
introduced by nonsymmetry is expected to be much less than the errors 
committed in past drag studies. The method presented here is far 
superior in terms of error than previous methods. 

Improvement of the errors in the present analysis could be made 
by a more accurate treatment of aerodynamics and more accurate 
measurement of accelerations and angles of orientation. The accuracy 
of angle measurements was about i 1° while the force measurements 
were accurate to i 5% for the data used. The aerodynamic description 
of the OVI-15 could be improved by employing a more accurate expression 
for the variation in Cp with I’j. The polynomial approximation employed 
in the analysis could be improved or the exact expressions could be 
employed at the expense of computer time. 

The results obtained for Cp as a function of angle of attack for the 
OVI-15 is of special interest because of the many drag analysis which 
have been performed on the satellite. For example, Champion, Marcos 
and Mclsaac, 1970; Marcos and Champion, 1972; Marcos, Champion, 
and Schweinfurth, 1971; have analyzed the accelerometer data of the 
OVI-15 to reveal a number of properties of the upper atmosphere, 
these analyses, accelerometer data was used only when the satellite 
was broadside into the flow. This instantaneous attitude would correspond 
exactly to the 0° or 180® spin axis orientation of the OVI-15. At this 
attitude Cj-j is equal to Cp as given in Figure 5 and would have a value 
of 2. 539^ for the case of infinite speed ratio and Ofj = . 64 and Pj = . 44. 

This value of Cj) is 15.4% greater than tlie value of 2.2 which was employ- 
ed in these analysis. Additional correction would have to be made if 
speed ratio effects were. taken into account. These corrections would 
result in substantial decrease in densities reported using a Cj-j of 2.2, 
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OVl-15 orbital data has been analyzed by Ching 1971 and King-Hele 
and Walker 1969. The King-Hele and Walker-analysis employed a 
constant of 2.2 independent of the satellite orientation. The reference 
area used by King-Hele and Walker was midway between the maximum 
of 2.578 nr^ and the minimum of nr^ shown in Figure 5. On the basis 
of the reference area employed by King-Hele and Walker, a Cq value of 
between 4. 543 and 3. 889 would correspond to the values of a j and Pj 
found in the analysis. 

The orbital decay analysis reported by Ching 1971 includes a factor 
to represent the changing aerodynamic drag properties of the OVI-15. 

The factor is based on the changing cross sectional area seen by the flow. 
The drag coefficient is considered constant while the reference area used 
in the analysis is changed by as much as 25%. A 25% correction factor 
is too large in view of the results of Figure 5 which show that the maximum 
change in Cj) A would be 16%. 

In addition, the effect these results have on past analysis of OVI-15 
data in particular, the results indicate that the assumption of C£) used 
in most drag studies have been too low. The value of Cj) = 2.0 or 2. 2 
which has been used for most past drag analysis is lower than could be 
expected for most shapes with Of j = .64 and Pj = .44. A sphere for 
example would have a Cjj , = 2.352 which is 7% higher than the 2 .2 
value often used. It should oe noted, however, that the results obtained 
here are for one satellite surface and one atmospheric composition. It 
is expected that other surfaces and other compositions should change the 
values of aj Pj result in a change in aerodynamic properties. 

More data must be collected before a firm value of oj and Pj can be 
assigned to a given gas and surface combination. Future work should 
be directed towards this goal. 
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TABLE I 


Orbit Number = 890 


Time of Peak F = 672 


I 

F(i, 1) 

GAMA(i, 1) 

F(i.2) 

GAMA(i, 2) 

TIME(i, 1) 

TIME(i,2) 

1 

317. 5 

137. 75 

320. 0 

134. 8 

647 

697 

2 

307. 0 

139. 50 

306.0 

133.25 

622 

722 

3 

291. 5 

141. 75 

291. 0 

132.25 

597 

747 

4 

270. 0 

142. 50 

268. 0 

130. 0 

572 

772 

5 

251. 0 

144. 0 

244. 0 

129. 0 

547 

797 

6 

225. 0 

145. 75 

214. 0 

127. 5 

522 

822 



Orbit Number = 893 


Time of Peak F = 

681 


I 

F(i. 1) 

GAMA{i, 1) 

F(i,2) 

GAMA(i, 2) 

TIME(i 

1) 

TIME(i, 2) 

1 

292. 0 

141. 75 

292. 0 

137.25 

650 


712 

2 

282. 0 

142. 75 

283. 5 

135. 0 

625 


737 

3 

264. 0 

144. 50 

263.0 

132. 75 

600 


762 

4 

243. 0 

146. 20 

239. 5 

132.0 

575 


787 

5 

218. 0 

147. 50 

215. 0 

130. 0 

550 


812 

6 

199. 0 

148. 20 

190. 5 

128. 0 

525 


837 



Orbit Number = 896 


Time of Peak F = 

675 



F(i. 1) 

GAMA(i, 1) 

F(i,2) 

GAMA(i, 2) 

TIME(i 

,1) 

Time(i, 2) 

1 

306. 0 

133. 5 

302. 5 

129. 5 

650 


700 

2 

302. 0 

13S.0 

291. 0 

127. 7 

625 


725 

3 

289. 0 

136.25 

275. 0 

126. 75 

600 


750 

4 

270. 5 

137. 5 

256.0 

125. 0 

575 


775 

5 

250. 0 

138. 0 

231.0 

123.0 

550 


800 

6 

226. 5 

141. 0 

200. 5 

120. 5 

525 


825 
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Figui'e 1. The gas surface interaction and the forces of the interaction. 
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Figure 3. Orientation of OVI-15 with respect to flow velocity U. 



Figure 4. Accelerometer output and angle of orientation gamma for 
orbit number 893. (From Fess and Young, 1969). 
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Figure 5. Constant values of the sum-of-the-squares -of-DEL^ over 

an acceptable range off '• and P- for data from orbits number 
890, 893, and 896. ^ ^ 





CHAPTER III 


AERODYNAMIC LIFT EFFECT ON SATELLITE ORBITS 

A paper submitted to the AIAA 13th Aerospace Sciences by 
G. R. Karr, J. G. Cleland and L. L. DeVries. Presentation was made 
during meeting held at Pasadena, California, January 20-22, 1975. 
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ASRODYMAKIC LVT VFBCT OH IATIU.ITI ORBITS^ 
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Unlvaralty of Alabama In Huntavllla 
HunCavllla, Alabama 35807 
and 

Dr. L. L. DaVrlai*** 

NASA- Marshal 1 Spaca Plight Cantar 
Huntsvllla, Alabama 


A bstract 

Numerical quadrature Is employed to obtain 
orbit perturbation results from the general per- 
turbation equatlona. Both aerodynamic lift and 
drag forces are Included In the analysis of Che 
satelllta orbit. An exponential atmosphere with 
and without atmospheric rotation Is used. A com- 
parison Is made of the perturbations which are 
caused by atmospheric rotation with those caused 
by satelllta aerodynamic effects. Results Indi- 
cate that aerodynamic lift effects on Che semi- 
major axis and orbit Inclination can ba of Che 
same order as the effects of atmosphere rotation 
depending upon Che orientation of the lift vector. 
The results reveal the Importance of Including 
aerodynamic lift effects In orbit psrturbatlon 
analysis. 


nonsynaetrlc pattern with respect to the velocity 
vector. Consider, for example, a cylinder or 
cone with an axis of symsetry along the dlracClon 
of the unit vector a. If such an object ware to 
travel through Che qtmosphara with a In the 
same direction as V, then only drag forces would 
result. However, If a were at an angle 8, with 
respect to then a lift force will arise given 
by 


-T- p C, A . g * 8 j X V 
^ sin e. 


( 2 ) 


where_ C, Is the lift coefficient of Cha object 
and A is assixaed the same as the X used In 
Che drag equation (Bq. 1). 


1. Introduction 

The perturbations of a satelllta orbit 
caused by Che Interaction of the satelllta with 
the earth's aCmosphara has bean a topic of con- 
siderable Intareat since orbital flight was pro- 
posed. The forces acting on a satelllta during 
Its passage through the atmosphere at speads of 
near 8000 m/sec are predominantly the drag force 
which we will define Co be Chat forca which acts 
parallel Co Che velocity vector, V, of Che sat- 
ellite with respect to the atmosphere. For a 
satellite having a drag coefflclant C_ and a ref- 
erence area A, Che aerodynamic drag force, 6, Is 
given by 

^D " “T “ ^ S ^ M 


where p Is Che density of the atmosphere. An- 
other aerodynamic forca which may act on tha 
satellite Is the aerodynamic lift force which 
we define as Che force perpendicular to velocity 
vector of Che satellite with respect to Che at- 
mosphere. Aerodynamic lift forces arise when a 
nonspherlcal aatelllte travels through the at- 
mosphere at an attitude such that atmospheric 
molecules are deflected by the satellite In a 
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The drag and lift forces given In equations 

(1) and (2) are both defined using the velocity 
with respect to Che atmosphere as distinguished 
from the Inertial velocity of the satellite v. 

The Inertial velocity Is defined here as Che ve- 
locity given by the orbital elements of Che oscu- 
lating orbit at Cha satellite. 

/u (1+ecosE) 

“V • h”- e co. hV 

where a and e are the osculating alemenCs of Che 
orbit and E Is Cha eccantrlc anomoly. Since Che 
atmosphere at orbital altltudas may have motions 
with respect to Inertial spaca, Cha Inertial ve- 
locity V Is not necessarily Ch« same as Che ve- 
locity of the satellite with respect Co Che at- 
mosphere 9. This difference Is due Co Che mo- 
tion of the atmosphere, ^ , at the satellite. 

The velocity ^ Is Chen moSlflad by Che atmospher- 
ic motion to give 

^ ( 4 ) 

Orbital perturbations undar Cha action of 
pure drag (l.a., 9, ■ 0)i with and without atmos- 
pheric motions wltn respect to Inertial refarenca, 
have received considerable attention for the pur- 
poses of (1) predicting Che orbit of a «atalllCa 
Into Che future (dee, for example. Ref. 1) and 

(2) deducing atmospheric properties from observed 
perturbations (see, for example. Ref. A). The 
Influence of nonzero lift (f. 4 0) has received 
little attention due to a number of factors which 
have been cited to reduce lift effects to negli- 
gible values. For example, random tumbling of a 
nonspherlcal satalllte Is oftan cited as a reason 
for neglecting lift since the lift vector would 
be randomly oriented and cand to average to zero 





undar such condttloM, Sacond, lift forcas ara 
thought to ba nrgllglbla for ■atallltaa baaad on 
tha argumant that tha aatallita la a poor raflac- 
tor of aK>laculaa alnca tha gaa aurfaca Intaractlon 
la asauned to ba Inalaatlc. Howevar, tha aaaiaap- 
tlon of Inalaatlc raflactlon la not varlfiad and 
tha poaalblllty atlll axlat* of algnlflcant lift 
coafflclanta. Whlla many of tha aarllar aatal> 
lltaa wera nearly apharlcally aynmetrlc and ran- 
domly tumbling (uncontrollad cltltudaa), many 
recant aatellltaa ara now apharlcal and may have 
large wlng-lika conf Iguratlona (auch aa Skylab) 
and generally require altitude control by either 
active or paaalva maana. Such aatellltaa may 
then experience lift forcaa of greater magnltudaa 
than experienced by peat aatellltaa. 

Another raaaon that lift forcaa have re- 
ceived llttla attention la tha poaalblllty that 
orbital parturbatlona that have bean cauaed by 
lift forcaa may have been wrongly attributed to 
other effecta, auch aa atmoapherlc motlona. At- 
moapherlc motlona that are not collnear with the 
Inertial velocity con give rlae to pure drag 
forcea with componenta perpendicular to the In- 
ertial velocity vector. One of the major motlona 
of the upper atmoaphere often Included In drag 
atudlea la that movement which la correlated 
with the rotation of the earth. The atmoaphere 
la as8«BDed to rotate In Inertial apace with 
about the aome angular velocity of tha earth. A 
pure drag aatelllte In an Inclined orbit would 
then experience componenta of the drag force 
which are perpendicular to the orbital plane. 
Theae forcea cauae the orbit plane to rotate In 
apace. The obaervatlon of ratea of change of 
the orbital plane of aelected aatellltea haa 
been uaed aa evidence of the rate of rotation 
of the upper atmoaphere under the aaaumptlon of 
pure drag only (aee Ref. 3). We will ahew that 
aerodynamic lift forcea can give rlae to orbital 
plane changea of magnltudea comparable to thoae 
cauaed by rotation of the upper atmoaphere, 

Thla la not meant to imply that there la no at- 
moapherlc rotation becauae there la atrong evi- 
dence to auggaat a co-rotatlng atmoaphere for 
the earth. We will ahow, hwever, that lift 
forcea of eeemlngly email magnitude can at leant 
produoe orbital plane changea that are of the 
some magnitude aa thoae aaaumed to be cauaed 
by atmoapherlc rotation. On thin baela we feel 
that lift effecta ahould be examined carefully 
In the rr'Mctlon of orbit perturbation data for 
the purpoee of deducing upper atmoapherlc mo- 
tlona. 


eccentrlcltlee. High eccentricity orblte eleo 
have the effect of concentrating all the aero- 
dynamic effecta at the perigee region alnce the 
atmoapherlc denalty dropa off exponentially away 
from perigee. 


II. Aarodynaaplc Lift 

In the free molecule environment at orbital 
altltudaa, aerodynamic lift and drag forcea are a 
direct function of the Interaction of tha atmoa- 
pherlc gaa moleculea with the expoaed aurfecea of 
the aatelllte. The characterletlca of the gaa 
aurfoce Interaction fer colllelon of upper atmoa- 
pherlc moleculea with aatallita materlela at ve- 
locltlaa of the order ;6f SOOO m/aec la not well 
undaratood. The gea i^v-faea Interaction lo de- 
acrlbed In a^M*^*'* manner ualng a modal devel- 
oped by Kar^ rn'whlch tbi« force acting on cn ele- 
mant of aurface axpoae-: to the flow la divided 
Into two componenta. With reference to Figure 1, 


< iwtacma 
Floe 



FIGURE 1. Illuatratlon of Forcea Acting Due 
to Gaa Surface Interaction 


the force acting on the aurfoce will have caeq>o- 
nenta F^ which la the force aaaoclated with the 
momentum carrlad by the Incoming moleculea and ‘ 

F* which la the reaction force aaaoclated with 
mSleculea leaving the aurface. The group velocity 
of moleculea leaving the aurface, D., la aaaumed 
to be aome fraction of the group velocity of the 
Incoming moleculea, D, due to poaalble Inalaatlc 
colllalona with the aurface. The relatlonahlp 
between and U la given by 


The purpoae of thla paper la to exmnlne the 
eeaentlal characterletlca of orbit perturbatlona 
that reault when a aatelllte haa lifting prop- 
ertiea. Comparlaona will be made with pure drag 
caeea where appropriate and comparlaon will alao 
be made with the effect of the atmoaphere ro- 
tating at the earth'a rotation rate. The baalc 
aerodynamic relatlona will be preaented and em- 
ployed In the general perturbation equatlone. 
Lift forcea will be divided Into two claaaea: 

(1) lift forcea In the orbit plane and (2) lift 
forcea perpendicular to the orbit plane. Numer- 
ical Integration of the perturbation equatlona 
waa made ualixg an exponential atmoaphere for pur- 
poaea of Illuatratlon of the aerodynamic lift 
effecta. Special emphaaia la given to the de- 
pendence of lift effecta on high eccentricltlea 
In order to expand on prevloua work valid at low 








Uj - ^ 1 - aj’ 0 (5) 

where th e proportionality term la expreaaed aa 
Vl ■ a.' In order that the parameter a. reaamble 
what la cuatomarlly termed the thermal accoomo- 
datlpn coefficient, A aecond parameter la lntr<^ 
duced to provide the direction of the reflected 
force. The angle of reflection law choaen lo a 
linear relatlonahlp between the angle of Inci- 
dence, 6, and angle of reflection, 6j, given by 

0j - I Pj ♦ (1 - Pj) e (6) 

where la the adjuatable parameter with P^ ■ 0 
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giving apecular reflection (6. « 0) and Pal 
giving diffusive reflection (0 an/2> Independent 

of 6). The mass flux Impinging on an element of 
aurface Is given by 

Ba-p0.ndA (7) 


where It Is the unit outward normal of the sur- 
face. The exchange of monentun that takes place 
at the surface gives the net force acting on the 
surface 


P « - (0 - t?j) p If . It d A (8) 

The total drag and lift forces acting on an ob- 
ject are obtained by Integrating equation 8 over 
the exposed surface of the satellite. Detailed 
results for spheres, cylinders, cones, and flat 
plates are given In reference 4 and 5. 


uncertainties In upper atmospheric density, accu- 
rate deternilnatlona of the gas surface Interaction 
from drag studies aloM.ls not feasible. Recent 
work by Reiter and Moe, ^however, on the analysis 
of drag and torque acting on a paddlawheel satel- 
lite provides the best values of gas surface In- 
teraction parsneters since atmospheric density Is 
eliminated In the analysis. Reiter and Moe, em- 
ploying a number of gas surface Interaction modela, 
concluded that the energy accoeiBodatlon a. Is 
In the range of .75 - .95. This would put>/l-a ' 

In the range of .5 to .22. The accoasodatlon co- 
efficients coming from their studies are high. 
Indicati ng Inelastic collisions. The value of 
yjl - O. Is seen to be significant even at high 
occoosnMatlsn coefficient values (for example, 
a. “ .99 gives >/l - a, ■ .1) which Implies that 
lift forces nay also Have significant values. 

The ratio of lift forces to drag forces act- 
ing on an object Is a nondlmenslonal measure of 
the magnitude of lift effects. For the flat plate, 
this ratio Is obtained from equations 8 and 9 


The above equations show that forces parpen 
dlcular to the velocity vector can only arise 
through the term In equation 8 . Since the 
magnitude and-^ direct Ion of P are functions of 
a. and P , the lift force acting on any object 
dlrectljrrelated to the gas surface Interaction 
parameters. This observation provides added mo- 
tivation for the study of llft-lnduced orbit per- 
turbations, since measurements of such perturba- 
tions could yield Information on the g«u-surface 
Interaction. 

For purpose of Illustration of lift- Induced 
orbit perturbations, flat plate drag and lift 
properties will be used, given by 

C_ B 2 sin 6 


271 - ttj sin 0^ Pj+(2-Pj) 0^ 

(9) 




.y l-a^' sin 

V 

. (! - f ) 9.1 


is 


1 1 -v'l-Oj cos 

TT P + ^2-P J 
2 J J 



s. - 


- 2jl - Oj' sin 0^ 

( 10 ) 


( 11 ) 


The lift to drag ratio given In the above equa- 
tion Is tabulated In Table I for four sets of gas 
surface Interaction parameters. The values are 
tabulated as a function of angle of attack and 
shew that (1) small angles of attack of a flat 
plate yield the highest L/D values, (2) specular 
reflection causes the highest L/D values (L/D - 
tan 0 , for specular reflection) and (3) highly 
Ineloltlc collisions (a. > .99, P » 1) cause 
lift forces which are almost 101 of the drag 
force. The values of a, • .75, P. ■ 1 were cho- 
sen to simulate the gas'^surface Interaction re- 
sults obtained by Reiter and Moe. This set of 
parameters produced L/D values of near 0.5 and 
this value of L/D will be used to Illustrate the 
orbital perturbation caused by lift forces. The 
Intermediate case (a. - .5, P. - .5) was chosen 
only to Illustrate the effect‘d of reflections 
other than specular and diffuse. T*'e L/D values 
for the Intermediate case are seen to be near 
unity at the small angles of attack. 


where 0 is the angle of attack tf the flat 
plate. For 0 >0 the plate Is edge-on Into the 

flow vrhlle fof 0^> tt/2, the plate has the out- 
ward surface normal In to the f low. Equation 10 
shows that the factor >/l - a.' plays an Important 
role In determining the lift^magnltuda. If the 
molecules have on inelastic collision with the 
surface, the velocity of reflection may be ex- 
tremely low, the value of a. would be near unity, 
and the lift very small. The perfectly elastic 
or perfect specular reflection (P ■ 0, a. - 0) 
would provide the highest possible values-^ of 
lift. 


Laboratory experiments at velocities corre- 
sponding to satellite velocities are not con- 
clusive on the gas surface Interaction to expect 
at orbital altitudes. Work by Hulpke^ for ex- 
ample tends to show elastic collisions while 
other experimenters have found diffusive type 
reflection (Ref. 9) or both (Ref. 6). Due to 


A typical satellite would not have the high 
L/D values as experienced by the flat plate which 
is an Ideal lifting body at orbital altitudes. 

For comparison. Table II gives the L/D values 
obtained from results given by Sentman (Ref. 10) 
for a typical satellite shape and typical gas 
aurface Interaction values. The shape Is a cy- 
linder with a conical end and having a total 
length of 4 times the diameter. The values of 
L/D are generally lower than the flat plate val- 
ues because this object maintains a high drag 
profile at all angles of attack. The peak value 
of L/D of .044 Is seen to occur at near 25*^ angle 
of attack. While this L/D Is much smaller than 
flat plate values, the lift force will be nearly 
5X of the drag force for this typical satellite 
« shape . 

Having nw established the likely existence 
of aerodynamic lift forces of magnitudes from a 
few percent of the drag force and greater, the 
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orbital parturbatlona that are characCeristtc of 
theta forcaa will be preaantad. In view of the 
fact that the aarodynaml: drag force la t)^ prin- 
ciple source of orbital perturbations, a lift 
force of even a few percent of the drag force la 
expected to be significant. In order to enhance 
the lift effects and to reduce cooputer tljsa, an 
L/D value of .5 was chosen for moat of the stud- 
ies. The caaq>uter sleulatlon employed a light 
aatelllta weight also for the purpose of enhanc- 
ing the aerodynamic affects. For these reasons, 
the results obtained are not likely to simulate 
any known satellite but are presented only to 
Illustrate the character of lift- Induced orbit 
perturbations. 


Table 1 


da 2 f _ .... 

7t ■ — 175 rnr/T * <“> 

(pa) (1 - a*) 



1/2 

1 

“T 


2g-e*) cos g 
1 - o cos B 


T 




1/2 


sin g 


N 


(13) 


Flat Plata Lift to Drag Ratios 


9s 

Specular 

Diffuse 
a. -.75, 

i'* 

Diffuse 
a, - .99 
‘ 

Inter, 
a, -.5, 

0 


6.3 

0.1 

1.0 

2“ 

28.636 

0.491 

0.100 

0.997 

5® 

11.430 

0.477 

0.099 

0.985 

K 

3. 671 

0.433 

0.097 

0.947 


3.732 

0.428 

0.094 

0.896 

25® 

2.143 

0.374 

0.087 

0.772 


1.428 

0.318 

0.078 

0.642 


1.000 

0.261 

0.066 

0.314 

K 

0.700 

0.203 

0.033 

0.392 

65® 

0.466 

0.143 

0.039 

0.276 

83® 

0.088 

0,029 

a.QM 

0.054 


Table II 

L/D for Typical Satalllte Shape 
(from Sentman, Ref. 10) 

Aaale of Attack L/D 


0 

0.0 

10 

0.037 

13 

C.041 

23 

0,044 

35 

0.033 

45 

0.023 

35 

0.013 


III. Orbit Perturbation gquatlons 


dl 

It 


r cos ( 8 

W 

(U •) (I 


± 



w 


(14) 


dQ 

dt 


sin 

J7? 


( Q Ml.} 
2. l/2. 


(U a)*'* (l-a‘) *'‘sln 1 


(15) 


dw 1 / s 2 sin g _ 

dt " oJT \ u / 1 - e cos g 


( 16 ) 



4 COS g) N 


(17) 


r sin (6 -t- uu) U (18) 

i j r - \ {/i 

(p a) (1-e^) tan 1 


where 

r 2 ,1/2 -1/2 

f m Kl-a^) (1 4 • COS g)J (1-e cos g) 

(19) 


The angle 8, the tnia anomaly, can be axpressad In 
terms of g, tha accentrlc anomaly, by 


The approach will be similar tc that taken 
by others In that the expressions for perturbing 
force will be substituted Into the general per- 
turbation equations for ths tins derivatives of 
tha osculating orbital elements a, a, 1, is, and 
0; tha semi-major axis, eccentricity. Inclination, 
argtssa nt of perigee, and right aseenslon of tha 
ascending node, respectively. Numerous good re- 
farencas provide details of tha derivation of 
these equations with specific application to drag 
forces (sea, for exmipla, Raf. 11). Tha forms of 
tha equations most useful for this work are given 
as 


of & page is 
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0 _ cos*^ ( f ^ (20) 

\ 1 - a cos g/ ' 

4tnd the radius, r, can also be expressed as a 
function of g by 


r ■ a (1 - a cos g) (21) 


Tha T cosiponant of force la tangent to tha orbit 
path and positive in tha direction of notion; the 
R component Is positive In tha direction of tha 



outward nonal to tha oxblt path and is In tha 
otblt plana; tha W coaponant la tha forca par> 
pandleular to tha orbit plana with positlva balng 
in tha same dlraction as tha orbit angular ao- 
mantisB vactor, Tha aarodynaaic fores acting on 
tha satallits will ba axprsssad in coaponants a- 
long M, T, and U. 

Considar a flat plate in orbit with tha 
orientation of tha plats held constant with re- 
ap* ct to tha M, T, W coordinate systaa. Tha unit 
vactor 1 will ba taken as tha normal to tha flat 
plate with componants in tha M, T, W diractiona 
glvan by 


•t 

■ cos 

•w 



m COB 

•w 

Sin 0jj^ 

■w 

- Bln 




Equations 1 and 2 give tha lift and drag forca as 
a function of tha ralatlva velocity V. Tha ra- 
lativa velocity Tf is first fouud with raspact to 
the R, S, W coordinate systan in which tha R di- 
rection is radially outward from tha gaocentar 
and S is perpend iculur to R and W with positlva 
being in the direction ox Tuition, An atmosphere 
rotating at tha earth's rotation rate, , is to 
b* considered the only atmospheric motloS for 
purposes of this study. The velocity fiald pro- 
duced by the atmospheric rotation is glvan by 


a -p ^eos sin 0^ (1 a eoa 6) 

* eoa 0^ cos 0^^ a sin sj 

(26) 

■ -p ^eos 0y eoa 0^ (1 a a eoa 0) 

- eoa 0^ sin 0||^ a sin oj 

\ - aln0,, 

Tha angls 6 callad for in tha drag and lift re- 
latlonshlps’can ba obtained in terms of tha orbit 
paraaetar and tha angles 0^ and 0^ from 

«. • | t | («> 

whare 

V ■ » i + (r 0-0^ r cos 1)8 +[0^ r eos(04w) 

sin i w] (20) 


. B. . r 


r n [cos i S - sin i cos(0-«w)w] (23) 


A A 

where S, W, and r are unit vectors, Tha inertial 
velocity of tha satallite is given by 




• A 

r R + 


r 6 8 


(24) 


The relative velocity Tf is »iven by aquation 4, 
Equations 4, 23, and 24 are substituted into 
equations 1 and ^ to obtain the components of lift 
and drag In the R, S, W coordinate system. The 
unit vector t will have componants a_, act and 
a^ In the R, S, W system. The transfonia^on 
from the N, T, W system to the R, 8, W system Is 
obtained from tha following relationship between 
the unit vectors 

H ■ — jr- [(1 + e cos 0) R - e sin 0 8 J 


A 

T 


(25) 


A 

sin 0 R 


e (1 -f e cos 


A ^ 

0) s] 


It is convenient to write the velocity ^ as a 
function of B and B, To do this, tha following 
relationships are employed 


s e 

r ■ a a sin B E 

, 1/2 . 

r 0 - a (1 - a^) E 

5/2 

r ■ (1 - a cos B)^ B 


(29) 


a^'^(l - a cos B) 


Substitution of these relationships into the fpree 
equations, and taking componants in tha R, S, T 
directions, we obtain tha following components of 
force R, S, and T. 

, _ A 

' • » <*-•«<-«> 

- ^ E CM e, ♦ ip } (») 


Using aquations 25 and 22, we find 
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(-S-) OS, »t e, ♦ v} 


(31) 
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(1 - • cos E) cot (9 4 ii)tln 1 


(»L ®®‘ ®a ♦ V I 


(32) 


C . A 

®D ■ ~"Tir~ • *L “ *TT" 




of aatalllte and n 


- V777. 


Thla concludaa tha development of the orbit 
perturbation relatlonahlpa. Equatlona 30, 31, 
and 32 along with the tr ana formation relatlona 
given In equation 25, provide the neceaaary re- 
latlona to write Che orbit perturbation equatlona 
In tha fora 


dE “ ^a^S.* ^D' ^W* ^NT* *• *» 

■Jg ■ *e^®L* ^NT' ** ** 

etc. (33) 


The atap alae on AB waa taken aa aaell aa 0,1 de- 
gree depending on the rate of variation of the or- 
bital elementa. The accuracy of the raau"ta la 
not expected to be good near the final atagea of 
decay alnca the orbital alaawnca are changing rap- 
idly. Brrora alao build up at 1 m eccentric It lea 
due to the high rate of change In d) which Invall- 
datea tha aaaunpClon of conatent orbital angular 
noeantum Implied In equatlona 29. 

The Integration of Che equ«.tlona atepwlae 
In B la a departure from prevloua work by Oook 
(Bef. 14) In which a aerlea expanalon waa eaiploy- 
ed to facilitate Che Integration analytically 
from 0 Co 2rr. We did not Cake thla approach for 
two reaaona: Fleet, Che reaulta obtained by Cook 

are only valid at low eccantrlcltlea (leaa than 
,2) and expanalona for high eccantrlcltlea were 
not found. Second, we were Intereated In the lift 
effecta on orbit decay which are not apparent If 
one Cakea, aa Cook did, equation 12 to Imply that 
lift forcea have no effect on decay of the aeml- 
najor axla, Bquatlon 12 ahowa that Che rate of 
change In a la, to flrat order, a function of 
only forcea tangent Co Che orbital path. For a 
non-rotating atmoaphere, drag forcea only would 
produce tangential componanta and lift forcea 
would not effect a. We find chat In a non- ro- 
tating atmoaphere, 0 ■ 0, lift forcea do affect 

the value of a over that of drag alone. Thla 
la clearly a aecond order effect which we find la 
comparable to the effect of atmoapherlc rotation 
on a for L/D ■ .5. Reaulta auch aa thaae are 
found by atepwlae Integration of Che orbit equa- 
tlona and would be loat If Integration from B ■ 0 
to 2n were done aaaumlng all parametera remain 
conatant over the Interval. 


V. Reaulta 


The reaulta are grouped Into three claaaea 
of lift orientation: (A) Lift forcea perpendicu- 

lar to the orbit plane, (B) Lift forcea In the 
orbit plane and In Che aame direction around the 
orbit, and (C) Lift forcea In the orbit plane 
which change direction at perigee and apogee. 

The flrat and laat type of lift orientation aervea 
to approximate caaaa In which aatellltea maintain 
an orientation with reapect to Inertial apace. 

The aecond type of lift orientation aervea to ap- 
proxlmace the earth oriented object. 


The atmoaphere denalty called for In t^a pertur- 
bation equatlona la given by 

P (r) - Pp e " 0'-hp)/H ^ 3 ^, 

where p la the denalty at perigee, h la Che 
height Bf perigee, H la Che acale helgRt at 
perigee, and h la the height of the satellite 
given by h - r - R_ where R^ la the radlua of 
the earth. ® 


IV. Hethod of Calculatlona 

The orbit perturbation equatlona were Inte- 
grated numerically ualng B aa the Independent 
variable, A modified Runge-Kutta fourth order 
process 'due to Olll^am discussed In ref, 13 was used. 
Double precision was used to maintain accuracy. 


A. Lift Perpendicular to the Orbit Plane 

For a flat plate with angles of orientation 
with reapect to the N, T, W coordinate system 
given by (0 • conatant ocher Chan aero, >0), 

lift forces will be produced In Che W direction 
which Is perpendicular to Che orbit plane. The 
effect on Che orbit of aerodynamic forces perpen- 
dicular to the orbit plane have been studied by 
Cook and Plljnner (Ref. 15) and Cook (Ref. 16) for 
the case of a satellite having pure drag with Che 
perpendicular forces arising due to the rotation 
of the atmosphere. For comparison, we took the 
cose of a non- rotating atmosphere with on aero- 
dynamic lift force of L/D - .5 and .1 acting per- 
pendicular Co the orbit. 

The orbit parameters most sensitive to forcea 
perpendicular Co the orbit plane la the orbital 
Inclination 1 and right aacension of the ascend- 
ing node n Aa pointed out In reference 15, Che 
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perturbations In 0 caused by aerodynsnlc forces 
perpendicular to the orbit are less than .IX of 
those caused by the oblataness of the earOh. The 
inclination i is not effected by earth oblate- 
noss to such an extent and Is therefore a sensi- 
tive indicator of aerodynamic forces perpendicu- 
lar to the plane. 

The characteristics of orbitnl inclinatloa 
changes caused by aerodynamic lift were found to 
be not much different from those caused by atmo- 
spheric rotation acting on a pure drag satellite. 
One difference is, of course, that the magnitude 
of caag>onent of drag force in the W direction 
caused by atmospheric rotation is propoi.tional to 
sin i whereas aerodynamic lift forces are inde- 
pendent of i. Therefore, inclination changes of 
near zero Inclination orbits must be caused by 
aerodynamic lift effects only, since atmospheric 
rotation can have no influence at zero inclina- 
tion. For comparison purposes, the lift induced 
rate of change in inclination with no atmospheric 
rotation was normalized with respect to the in- 
clination rateg of change for a pure drtg satel- 
lite at i ■ 45 in an atmosphere rotating at the 
earth's rotation rate. These results are pre- 
sented in Flgura 2 for L/D » .5 md L/D - .1 
with Bp - ,2 m^/Kg and Bj^ » .1 m^/Kg. 



imiiai IcetnifKiiv. • 

FIGURE 2. Out-of-Plane Effects On Orbital 
Inclination 


The notation In used to signify that the lift 
force is in the W direction. Figure 2 shows that 
the rate of change in inclination for L/D • .1 is 
over twice that that would be caused by atmo- 
spheric rotation. The curve for L/D > .5 shows 
that the dffact is proportional to L/D. Figure 
2 also shows a comparison of the atmospheric ro- 
tation at 1 ■ 90^ with respect to that at 1 ■ 45^. 
The effect is proportional to the ratio of the 
sines of the angles as would be expected. 

Figure 2 shows an increase in the rate ratio 
as a function of a. The reason for this increase 
is due to the decrease in di/dt for the atmo- 
spheric rotation effect as e is increased since 
Che inertial velocity near perigee Increases with 


a. The atmospheric rotation effect on 1 de- 
creases since It Is proportional to the ratio of 
V to V which decreases with eccentricity. The 
alrodynsmlc lift force is not affected by e since 
L/D la independent of velocity. 

B. In-Plane-Lift Forces Constant in Direction 

Consider now a flat plate with angles of 
orientation such that ■ 0 and « constant 
(other than zero). This orientation would pro- 
duce a lift force that is in the N direction and 
constant in sign around the orbit. Two cases will 
be considered, (1) L/D - .5 with the lift force 
directed always in the positive N direction and 
(2) L/D - .5 with the lift force directed always 
In the negative N direction. As pointed out above, 
aquation 12 shows that these lift forces do not 
cause a change in a to first order. 

We Investigated the change in a and e (4 a/ 
t t and ^ e/A t) at a point near perigee, but al- 
ways after perigee, using the numerical methods 
described. The results are given in figures 3, 4, 
and 5 in which the effect of lift and the effect 
of atmospheric rotation are both compared to the 
affects caused by pure drag only with no atmospher- 
ic rotation. The notation L^ Is used to signify 
that the lift force is in the N direction. 

Figure 3 shows that 6 a/A t is positive or 
negative depending upon the orientation of the 
lift vector. The magnltuda of the effect is de- 
pendent on the eccentricity, being larger for 
smaller values of e. Also plotted in Figure 3 
is the ei'fect on A a/A t due to a rotating atmo- 
sphere ^i - 45 ; . The magnitudes of the two ef- 
fects are comparable and are seen to be of order 
5X of the drag- nonrotating- atmosphere effect a- 
lone. The lift effects are second order effects 
as pointed out above while the atmospheric ro- 
tation effects are first and second order. The 
first order effect due to atmospheric rotation 
comes in because the force in the T direction ^s 
decreased by the wind for i between 0 and 90 . 
Figure 4 shows that the effect on e has much 
the same character as the effects on a. 
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FIOURE 3. In-Plano-Llft Effects on Beti.l-Major 
Axis Decay Betas 
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FIDURE 4. In-Pl«n*-Ll£t Effect# on Eccentric- 
ity Decay Rates 


Figure 5 shows the dependence of A a/A t on 
perigee height for two eccentricity values e - .2 
and e - .8. Below a perigee altitude of 130 Km, 
the rotating atmosphere effect Is seen to tend to 
reverse while the lift effects become larger and 
do not change sign as does the atmospheric rotation 



FIGURE 5. In-Plane-llft Effects on Instanta- 
neous Semi- Major Axis Decay Rates 
as a Function of Perigee Height 
The effect of constant lift on the orbital 
lifetime Is shown In Figure 6, The lifetime Is 
computed from time of initial conditions until the 
altitude becomes circular. Although the computation 
technl<)ue was not deslmed to handle zero eccentric- 
ities, little error Is Involved since the orbit 
eccentricity usually would become zero iron Ini- 
tially high values only very near final decay. The 
lifetimes arc normalized with respect to the life- 
time of a pure drag satellite having the aame drag 
as the lifting satellite. The atmosphere was taken 
to have zero rotation In order to clearly show the 
effect of lift. The results show thst In-plsne lift 
affects tend to decrease the lifetime for high ec- 
centricity orbits, of the direction of 
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FIGURE 6. In-Plane-Llft Effects on Satellite 
Lifetimes 


the lift vector H (assumed to be constant In sign 
around the orbit). At low eccentricities, the ef- 
fect on lifetime Is less but there Is separation of 
the effect according to the sign and magnitude of 
the lift vector. B»e results at low eccentricities 
contain some error due to the rapid rotation of the 
orbit In space. For this reason, the magnitudes of 
the Influence on lifetime may be In error at low 
eccentricities by an unknown amount. The curve Is 
presented here In order to show the clear tendency 
of lift effects to reduce lifetime as eccentricity 
la Increased. 

With reference to equation 16, the effect of 
constant lift force around the orbit should have a 
significant effect on i since the N component of 
force Is multiplied by cos E. This was also noted 
and evaluated by Cook (Ref. 14) for low eccentrici- 
ties. Using the same aerodynamic values used by 
Cook, we evaluated J) for low and high eccentrici- 
ties and compared the results as seen In Figure 7, 

The comparison with the Cook result Is seen to be 
good at the lowest eccentricities used In our nuMr- 
Ical work. In addition, our results show that it 
for higher eccentricities continues to decrease with 
Increasing e. 

C. In-Plana-Llft Forces Which Change In Sign at 
Perigee and Apogee 

Consider a flat plata flcwn such that the 
angle of attack In the HT plane changes dlscontlnu- 
ously In sign, but remains constant In magnitude, at 
perigee and apogee. This special case Is of Interest 
since It has the effect of approximating a space- 
craft with attitude controlled with respect to In- 
ertial spase. TWO cases arise; (1) lift In positive 
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FICURE 7. Argument of Perigee Rete of Change 
with Lift and Drag 


direction of N before perigee and In negatl\e di- 
rection after perigee and (2) lift In negative di- 
rection of N before perigee and In positive di- 
rection after perigee. These types of lift his- 
tories have the effect of making the N sin E terra 
In the equation for de/dt (Eq. 13)elther always 
negative, case 1, or always positive, case 2. 

Case 1 Is the same as that considered by Cook In 
reference 14». Instantaneous values of H aM t 
and A a/A t at a point Just past perigee for both 
positive and negative values of lift are given In 
figures 3 and 4. Values of A e/T where T Is 
the orbit period were com; Jted and found to agree 
well with the results published by Cook and will 
not be presented here. Of special Interest, how- 
ever, was the effect of the discontinuous lift 
cases on orbit lifetime which was not considered 
by Cook. 

The lifetime of the orbit was found to be 
strongly affected by the discontinuous lift cases. 
The results are given In Table III In which the 
lifetime is normalized to the lifetime of a satel- 
lite having drag of the same magnitude as the lift- 
ing satellite and with no atmospheric rotation. 

The lifetimes for case 1 lift histories were seen 
to be Increased by 1.7 to 3.5 times over the drag 
only lifetimes. Only three lifetimes were coa- 
puted due to the long computing times required as 

Table III 

Discontinuous Lift Effect on Lifetime (L/Dc,5) 


(Lifetime) /(Pure Drag Lifetime) 


e 

Case 1 

Case 2 


N sin E is Neg. 

N sin E is Pos. 

.1 

1.7 

.52 

.2 

3.2 

.48 

.3 

3.5 

.41 


the eccentricity was increased. The case 2 lift 
history Is seen to reduce the lifetime to about 
half compared to the drag-only case. Careful anal- 
ysis of our results shows that the case 1 lift 


history has the affect of reducing the perigee de- 
cay rate of the orbit, thereby reducing the drag 
effect, whereas the case 2 lift history causes the 
perigee height to decrease faster than the pure 
drag case. This effect Is best understood by con- 
sidering the relationship for perigee radius, r , 
given by ^ 

r^ - (1 - e) a (35) 

The time derivative of r Is 
P 


dr 

_E 

dt 


(1-s) 


da 

dt 



(36) 


The values of da/dt and de/dt are negative for pure 
drag and for both case 1 and case 2 lift histories. 
The action of case 1 and case 2 lift histories Is 
to modify the decay of perigee height with respect 
to the pure drag case In the following manner: 

dr / dr \ 

■ (1 i *) ( ) + (higher order terms) 

(37) 

where X Is a number of magnitude .1. The plus 
sign is taken for case 1 and the negative sign Is 
taken for case 2. The Subscript o refers to the 
pure drag values of r . 


VI. Conclusions 


Satellite lift forces were shown to be of the 
order of a few percent of the drag force for an 
ordinary satellite and may range up to near unity 
for a satellite designed to have high lifting forces. 
Reference to figure 2 concerning the lift effects 
on the orbit Inclination leads us to conclude that 
an L/D of only .01 would produce an orbital in- 
clination rate of change of a magnitude nearly 25% 
of that attributed to atmospheric rotation. Clear- 
ly, then, an attitude stabilized satellite with a 
lift force consistently perpendicular to the orbit 
plane could produce large changes In the Inclina- 
tion which may be wrongly attributed to high ve- 
locity upper atmospheric winds. 

The results obtained concerning the effect 
of different lift histories on the lifetime of the 
orbit are conclusive on three points: (1) constant 

lift directed either positive or negative along N 
about the orbit causes a reduction in lifetime at 
high eccentricities, (2) a discontinuous lift with 
change in sign at perigee and apogee, such as to 
have positive lift before perigee, causes Increased 
lifetime over a drag-only satellite, (3) a discon- 
tinuous lift change of the other type. In which 
the lift Is negative before perigee, causes a re- 
duced lifetime of the satellite. These conclusions 
are apparent from the results given in figure 6 
and those given In Table III. 

Due to difficulties we had In the numerical 
procedure near final decay of the satellite, the 
numbers on lifetime may be in error. For this 
reason, we stress that the results merely provide 
the character of the orbital perturbations due to 
lift and may be in error In absolute value. We 
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found It difficult to put • realistic aatlaata of 
error on the reaults obtained. All wock waa car- 
ried out In double preclalon and tha atap alta waa 
dacrcaaed until no changea In pravloua reaulta were 
found. We did find variations from tha smooth 
curves of figures 3, 4, and 5 of the order of 1 to 
SX. The results are considered accurate enough to 
Indicate tha controllability of satellite orbits by 
lift. 

The results show the influence of satellite 
lift on the orbital elesMnts and also open the 
possibility of utilising satellite lift to non- 
propulslvely control tha satellite orbit plane and 
orbital lifetime. A flying spacecraft would have 
an advantage over the propulsive satellite If the 
aerodynamic shape could be Incorporated Into the 
design without Increased weight. The weight of tha 
propulsion system and Its fuel could then be used 
as payload. 

Future work should be done on the observation 
of lifting effects on existing satellite orbits In 
order to confirm the results obtained here. In 
addition, analytic work and further numerical 
studies need to be done on the influence of lift 
on the orbital elements. 
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CHAPTER IV 


DUAL FALLING SPHERE DETERMINATION OF DENSITY AND 
TRANSITION FLOW PARAMETER 
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G. R. Karr. Presentation was made during meeting held In Washington, 
D. C. , January 30-February 1, 1974. 
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Abstract 

A new approach to the analysis of falling 
sphere drag data Is described In which the data 
from two trajectories through the same region of 
the atmosphere are analyzed simultaneously. The 
analysis provides Important aerodynamic Infor- 
mation which Is used to obtain an Improved value 
of atmospheric density. The technique Is applied 
to a set of falling sphere data In which a sphere 
transition-flow parameter and atmospheric density 
results are obtained In the 80-120 km region from 
published data for falling spheres over Kwajaleln. 
Another set of data for a falling sphere test over 
Wallops Island is also analyzed with comparable 
results. 


Introduction 


The falling sphere technique has been the 
prim*, source of measurements of atmospheric 
density and temperature In the Important altitude 
range of 80 to 120 km. In particular, the results 
from three falling sphere experimental groups 
provided the Information used to supplement the 
1962 U. S. Standard Atmosphere.*'^ These groups 
were Peterson, Hansen, McWatters and Bonfantl,’ of 
the University of Michigan; Falre and Champion" of 
the Air Force Cambridge Research Laboratories; and 
Pearson* of Australian Weapons Research Establish- 
ment. The results obtained by these groups are 
summarized In the 1966 supplements to the U. S. 
Standard Atmosphere.^ The method of analysis In 
all these experiments was to first measure the 
acceleration, a, acting on the sphere from either 
trajectory analysis or accelerometer readings. The 
drag equation 

Drag - j 0 V*Cd A - ma (1) 

Is then employed to deduce the density, p, where V 
is the velocity, Cp Is the drag coefficient, A Is 
a reference area and m la the mass. The mass, 
acceleration, velocity, and area are all measured 
quantities, leaving only p and Cd as unknowns In 
equation (1). A table of Cp as a function of 
Reynolds number, R^, and Mach number, M, Is usually 
employed to obtain a value of Cp. A density value 
Is then determined by solving equation (1), Since 
Cp Is a function of Reynolds number and Mach 
number which are density and temperature dependent. 


*Thls work was supported by the National Aero- 
nautics and Space Administration under Contract 
NAS8-28248 through NASA/MSFC, Huntsville, 

Alabama. 
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the analysis of the drag data usually Involves a 
number of Interatlons until there Is convergence 
to a final density value. 

The work reported here is the result of an 
examination of the falling sphere technique for 
the purpose of proposing possible Improvements, 
particularly In the area of aerodynamics. As with 
all drag deduced density exper Iments, , the values 
employed for the drag coefficient in the data 
reduction is a major source of error. The possi- 
bility of error due to drag coefficient Is largest 
In the 80 to 120 km region, due to the passage of 
the sphere through various aerodynamic regimes. A 
typical falling sphere trajectory will be In a 
free molecule flow regime at high altitude and will 
pass through transition flow Into contlnum flow at 
low altitudes during descent. In addition, the 
speed of the falling sphere may pass from subsonic 
to supersonic and back to subsonic during a typical 
traj ectory. 

A considerable improvement In knowledge of 
sphere drag coefficients has been made recently 
through a comprehensive experimental program at 
ARO which was sponsored by AFCRL for specific 
application to the falling sphere program. The ARO 
work reported by Bailey and Hiatt* covers a 
velocity range from 0.1 to 6.0 In Mach number and 
a Reynolds number range from 20 to 100,000. While 
this data Is extremely useful, there Is still a 
lack of accurate Information for the near-free 
molecule drag coefficients which would correspond 
to Reynolds numbers below 100. 

In the work reported here, two new methods 
of falling sphere data analysis are proposed 
which should find application In the low density, 
high altitude region of the atmosphere where 
accurate aerodynamic Information Is still lacking. 
The proposed methods Involve the simultaneous 
analysis of the trajectory of two spheres 
travelling through the same region of the atmo- 
sphere with different velocity. As an Illus- 
tration, the falling sphere data of Peterson,* 
et al. and Peterson and McWatters* Is analyzed 
using one of the proposed techniques. 

B asic Theory of Dual Falling 
Sphere Experiments 

Consider an experiment In which two sphere 
drag measurements are performed In the same 
region of the atmosphere but at different 
velocities. 



»1 - i p V Cdx S 

(2) 

t>2 • i P V Cd2 A 2 

(3) 

Since the apherea are In the aame region of the 
atmoaphare, the valuea of denalty are the aame for 
equation (2) and (3). Aaaume alao that the 
dependence of Cq on p, V and a third quantity B, 
which will be dlacuaaed later, la well underatood. 
Then, In general 

Cdj (P. B, Vj) ^ Cpj (p, B, V 2 ) 

(4) 

and equatlona (2) and (3) repreaent a 1 
equatlona In the unknowna p and B. The 
p and B are found by aolvlni equatlona 
almultaneoualy. 

net of two 
valuea of 
(2) and (3) 

P " p (Oj, D 2 , V^, V 2 , Aj, A 2 ) 

(5) 

B - B (D^, D^, V^, V^, A^, A^) 

(6) 


In the proceaa, not only la a value of danaity 
obtained but alao the two drag coefflclenta are 
determined for the two caaea. 


The baalc theory of dual falling aphere 
analyala dependa on the fact that the drag 
coefficient la not Independent of the velocity and 
can be tn-ltten ae the function of two other 
paraatetera at moat. The dependence of the drag 
coefficient on denalty and the quantity B nuat be 
known In order to write equatlona (5) and (6). The 
quantity B will be aeen to be the temperature 
for the flrat caae dlacuaaed below and the 
tranaltlon flow parameter for the aecond caae. 

If the drag coefficient were a function of more 
than two unknowna, for example 

Co ■ Cp (p, By, B^, B^, • • • Bn) (7) 

then two poaalbllltlea exlat. One could conalder 
multiple falling aphere experlmenta dealgned ao 
that m > 1 drag meaaurenent %rould provide Infor- 
mation neceaaary to Invert the m -f 1 drag 
equatlona, giving 

p-0 (Di,D 2 ••• Dnfi.Vi, Vj. 

®1 ■ ®1 <®1* ••• ®m+l» Vi ) 


leaa than In alngle falling aphere analyala where 
all the B^ quantltlea nuat be aaaumed. The approach 
taken In the follo«rlng dlacuaalon la to reduce the 
unknowna In the drag coefficient to the denalty and 
one other quantity called B In the above. The 
potential of the dual falling aphere technique 
appeara to be greateat In the high altitude region 
In which the flow la free molecule and the drag 
coefficient la very aenaltlve to the atmoapherlc 
tamperature aa dlacuaaed In the next aectlon. 

Temperature and Penalty 
Determination In Free Molecule Flow 


In the free molecule flow regime the drag 
codfflclenta of a aphere la atrongly dependent 
upon the apeed ratio for apeed ratloa of order 
unity and leaa. The apeed ratio la diiilned aa the 
ratio of the aphere velocity to the thermal 
velocity of the gaa given by 

S - W/y/l RT (9) 


the free molecule aphere drag coefficient la 
obtained from free molecule theory to be 



( 10 ) 


where K la a factor of order unity dependent on 
the gaa aurface Interaction. 


The drag coefficient la found to be Inde- 
pendent of denalty In a free molecule flow. The 
nature of the dependence of on S la better 
llluetrated by the expanalon of equation (10) for 
the caaea of large and amall valuea of S. The 
reaulta are 


Cd 


fa 


ml [il + X S--2- s’l (1+K) 
IT [3 S 15 240 J 


( 11 ) 


S <. 1.25 




( 12 ) 


S > 1.25 


B 2 “ B2 (Dj, D2, ••• Vj ) 

a 

®m " B, (Dx, •••• D^i, V^, ) 


(8) Equation (11) ahowa that tenda to Infinity 

aa S tenda to zero while equation (12) ahowa that 
^Ofn tenda to 2 (1-fK) aa S tenda to Infinity. 
Valuea of S of order unity and leaa are aeen to 
cauae the greateat variation In 


The aecond poaalblllty, and the one more 
practical to conalder, would be to aaaume all but 
two of the quantltlea needed to determine the drag 
coefficient In equation (7). An error la of courae 
BMde depending upon the accuracy of the 
aeaumptlona but the total error la expected to be 
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since falling aphere velocltlea at high 
altltudea typically correapond to apeed ratloa 
of unity or leaa, a dual falling aphere analyala 
may be feaalble and fruitfully applied In thla 
region. Conalder two aphere drag meaauramenta 
at free molecule condltlona 



■l ■ 2 ** V (Sj, Kj) 

- I PVj* (T, Vj, Kj) (13) 

-2 *2/^2 “i PV/ Cd,, (S2. Kj) 

■i‘’V <1)f. '^2' '^2> 


where the speed retlo hes been written es e 
function of tenpereture which nust be the same for 
both neesurements. Since free molecule drag 
coefficient Is independent of density, the density 
cen be eliminated In the above equations giving 


*1 

(T. Vi. Kj) Ai 


(15) 


2''*%, (T, Vj, Kj) Aj 

Equation (15) now contains the unknowns T, K^, and 
K 2 > If the two spheres have the same surface 
properties and V^ls not much different from V 2 , one 
would expect the gas surface Interaction to be the 
same for both spheres. Therefore, 


however, la lndepen 4 ent of the density determi- 
nation and becomes more accurate at higher 
altitudes where free molecule conditions prevail. 

A survey of published falling sphere data did 
not produce data of the type needed for an example 
analysis of the free molecule type. The proper 
data could be obtained by launching two spheres at 
near the same time but with different velocities. 
Another technique would be to track the sphere 
both during ascent and descent since the veloci- 
ties would be different due to drag effects. Data 
of the latter type Is available but only at lower 
altitudes where the flow Is transition rather 
than free molecule. Analysis In the transition 
regime Is considered In the follotrlng section. 

Dual Felling Sphere Analysis 
In th'. Transition Realms 

As a falling sphere passes Into regions of 
greater density, the dreg coefficient nust change 
from a free molecule value (2 and greater) to a 
continuum value (1 and less). The flow regime 
between the limits of free molecule and continuum 
Is termed the transition flow regime. Since no 
theoretical expression Is available which can be 
accurately applied to the transition flow regime, 
empirical and seml-emplrlcal relationships are 
comsonly employed. One such relationship given by 
Matting* has application to drag coefficient 
determination In the near-free molecule side of the 
transition regime. The expression Is given as 


Kl - K2 - K 


(16) Cd - CDc + (Ccf, - CDc) e-E/Kn (19) 


and Che common factor (l-*-K) can be eliminated 
from equation (15), leaving a single equation In 
the unknown temperature T. Therefore, from 
equation (15) 

T " T (mj a, , t ,2 * 2 » Vj, V 2 , Aj^, A 2 ) (17) 

and either equation (13) or (14) may be used to 
obtain 

p - p j^m a, (T, V, K)] (18) 

Notice that the value of K becomes Important In 
determining the value of density but Is not 
required In determining temperature. 


where Is the continuum drag coefficient, Kn Is 
the Kundsen number, and E Is Che parameter which 
must be determined from experiment. The above 
equation Is seml-emplrlcal based on a first 
collision analysis of near-free molecule flow. 

For application to falling sphere analysis 
the quantity E/Kn can be expressed ss a function 
of density since 



where r Is the sphere radius. Therefore, write 

E/Kn H C 3 p r (20) 


Discussion of Proposed Free 
Molecule Analysis' 

The sbove procedure would provide a more 
accurate measurement of temperature at high 
altitudes than that provided by single sphere 
experiments. In single falling sphere experiments, 
the temperature la deduced from the density values 
by Integration of the hydrostatic equation 
beginning at the high altitudes. As discussed by 
Bartman, Chang, Jones and Liu,* this method of 
datarmlnlng temperature Is subject to large errors 
at Che high altitudes. The dual sphere analysis. 


where C 3 will be ten.ed the transition flow onset 
parameter. Since and Cj)^ are functions of 

velocity, tenperacure, and the gas surface Inter- 
action, the functional dependence of the 
transition flow drag coefficient Is written from 
equations (19) and (20) 

Cd - Co (P, C 3 , T, V, K) (21) 

Since Cq contains four unkno%m parameters, 
t%K> parameters must be assumed In order to perform 
dual luUy sphere analysis In the transition (falling) 
.uglme. The value of temperature will be asstimed 
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to b« given by the ueuel acthode of intcgretion of 
Che hydroetetic equation since Bartaan,* et al. 
point out that the teaperature deduced by the 
Integration aethods becomes aore accurate at the 
lower altitudes. The aecond parameter to be 
aesuaed In the analysis la Che gas surface Inter- 
action paraaeter, K. The value of K • 0 will be 
employed In Che analysla for reasons dlscuseed 
later. A gas surface Interaction In which 
molecules are reflected In Che specular direction 
corresponds to K • 0, Independent of the degree 
of accoamodatlon to the surface temperature. 
Therefore, C3, la chosen ne the unknown paraaeter 
which will be determined In the analysis In 
addition to Che density p. The value of C3 la not 
well establlahed due to a lack of experimental 
results In the near-free aiolecule regime. The 
analysis will then help eatabllsh the value of 
this Important parameter for use In future experl- 
aenta. 

Having chosen the unknown parameters, the 
method of analyala la alallar to the free molecule 
flow analysla. Consider two drag measureaents 
In the transition regime at i.iie same region of the 
atmosphere but at different velocities. 

Di - p Vi» [0 d^j+ (CDj^-Cd^I )e-C3Pf]A 
1 r 

Dj ■ 2 P V2 * [Cd^ 2'" ^‘^Dfa2- <=Dc2 > ^ 

where the spheres are taken to be the sane size so 
that A3 ■ A2 ■■ A and r^ ■ t2 ■ r. The quantity 
exp (-C3pr) can be eliminated In the above 
equation giving a single equation In the unknown 
density p. 

(Di/>?A) (Cd£,2-CDc2)-<D24v|a) (CDf.j-CDei) 

p- f (23) 

‘^Dcl %ml 

and C3 Is found from either of equations (22) once 
a value of p Is determined from equation (23). 


C3 



(D/| V* A) ^ - Cp^ 


(24) 


The transition flow snalysls has been 
successfully applied to five sets of falling sphere 
measurements over Kwajaleln made by the University 
of Michigan group In 1963 and 1964 as described In 
the next section. 

Analysis of Kwajaleln Palling 
Sphere Measurements 


The University of Michigan falling sphere 
aeasureaents consist of data taken both during 
aacent and during descent for one of the three 
spheres ejected from the rocket during the aecent 
phase. The spheres were made of Mylar Inflated 
with Isopentane having an Inflated diameter of 


0.66a and a aass of 50 grams. The velocity- 
altitude history of a typical flight Is shown In 
Figure 1. In this particular flight there la an 
overlap of ascent and descent dsta In Che region 
between 90 and 110 km. Data of the type shown In 
Figure 1 are suitable for analysis as a dual fall- 
ing sphere experiment using the method outlined 
above. 


Unfortunately, not all flights had regions of 
overlap. Of the 13 successful flights over 
Kwajaleln, only six have any overlap. The sounding 
number and the region of overlap of the ascent and 
descent are the following: 


Sounding 2; 
Sounding 3; 
Sounding 8; 
Sounding 12; 
Sounding 13; 
Sounding 14; 


100 CO 102 km 
99 to 102 km 
104 to 109 km 
99 to 104 km 
96 to 107 km 
at 102 km 


The plots of these data are given In the reference 
by Peterson,’ et al. The data used for the work 
reported here was obtained In Che fora of computer 
output of Che University of Michigan analysis, s 
sample of which Is also given In reference 3. The 
temperature data obtained In Chat analysis was 
used directly In this work after smoothing over Che 
data In the region from 80 to 120 km using a third 
order least squares polynomial. The drag term 
required In equations (23) and (24) was obtained 
using the published values of Cp and p In the 
following relationship 


- (p c,,) 

^ j “ published values (25) 

2 '' ^ 


Due to noise In the data, the application of 
equations (23) and (24) could not be made point 
by point. A smoothing of Che data was than 
performed over the region of 80 to 120 km employ- 
ing a second order least squares polynomial fit. 
This meant that the ascent and descent 
trajectories were smoothed over a distance of about 
20 km. 

The Cp, values used were obtained from either 
equations (ll) or (12) and Cp^ values were obtained 
uelng a polynomial fit to the data given In 
Reference 6 for the highest Reynolds numbers. The 
mean molecular weight needed to determine speed 
ratio or Mach number was obtained from s polynomial 
fit to the mean molecular weight given In the 1962 
U. S. Standard’ for the region between 80 and 120 
km. 


Discussion of Results for 66 am Spheres 

Values of were calculated In Che region 
of overlap for each cf the sound Inga. Due to the 
smoothing operations only one value of could 
be obtained from each sounding. Soundings number 
8, 13, and 14 were analyzed using all the data 



points provl4«<l in the 80 to 120 ka region. The 

reeulte ere 

Cj( 8 ) •» 3.8200 X 10* m*/kt 
CjdB) - 3.2352 X 10* ■*/kg 
Cj(14) - 5.2345 X 10* a*/kg 

The dete for soundings 2 end 3 were noisier 
then the rest end s vslue of C 3 could only be 
obtained after ellainstlng a nuaber of the aoat 
divergent points. The results were 

C3<2) - 5.852 X 10* m*/kg 

Cj(3) - 15.924 X 10 * eVkg 

The noisy data associated with these results is 
likely the cause for the values being higher than 
for Che other three. 

The data for sounding 12 was smooth yet a 
solution for C 3 could not be obtained. It should 
be noted, however, Chat sounding 12 was also cause 
for concern by the Michigan group because of the 
anomoulous behavior which they felt might be due 
to a small leak in the inflated sphere. 

Results From 7 In. Sphere 

Additional overlapping falling sphere data was 
obtained for an accelerometer instrumented 7 in. 
sphere experiment over Wallops Island flown in 1961.- 
The data is published in NASA-CR-29 by Peterson and 
MeUattera^ consisting of results obtained from the 
first tests of the accelerometer system. The 
method of analysis was the same as used on Che 
66 cm data. The result obtained was 

C 3 (7 in) - 1.5583 x 10* m^'kg 

This result is within a factor of two of that 
obtained for soundings 8 , 13, and 14. The factor 
of two difference may be due Co the different 
surface properties of the 7 in. as compared to the 
66 cm. Also, the 7 in. sphere enters the 
transition regime at a lower altitude than does 
the 66 cm sphere due to its smaller sire. The 
different molecular compoalclon at lower altitudes 
may cause a change in C 3 . These questions could 
be answered by analysis of more 7 in. trajectory 
data. 

Atmospheric Density Calculations and 
Discussion of Results 

After determining s value for C 3 , Che drag 
data taken in Che original experiment may be 
reanalyzed to obtain new values of density employ- 
ing equation (22). An Iterative technique was 
employed using the unsmoothed tamperature and 
drag values given in Reference 3. A comparison 
was made of the density values given in reference 
3 with the density values obtained using . _ 


C 3 - 4 X 10*m*/kg which represents an average of 
the results for soundings 8 , 13, and 14. A point 
by point comparison was made for aach of the 
aoundlngs 8 , 13, and 14 and the average was 
obtained. The reaults are presented in Figure 2 
which shows that Che densities calculated using 
the methods described In this paper produce 
generally higher values than obtained in the 
original experiments. 

Figure 3, obtained from reference 2, shows 
the mean values of density obtained in all 13 of 
the original experiments as compared to the 1962 
U. S. Standard.' This figure shows that the 
original analyals resulted in a nearly lOX lower 
density value above 100 km Chan gives in refer- 
ence 1. Application of the results of the current 
work shown in Figure 2 would cause the density to 
be nearly equal or somewhat greater than Che 
U. S. Standard above 100 km. Both methods of 
analysis give a higher density in the 90 to 100 km 
region than that given in the U. S. Standard, but, 
neither analysis is correct in this region as 
discussed in the following: 

In the 90 to 100 km region, the Knudsen number 
is of order .1 as is shown in Figure 2. Transition 
flow is considered to be within the limits of 10 
Co .1 which means that in the 90 to 100 km region 
the flow is near continuum rather Chan near free 
molecule. For this reason, equation (19), which 
is derived on Che basis of near free molecule 
theory, is likely not valid in this lower region. 
The results could be improved by employing a more 
valid relationship. The method of analysis would 
remain the same however. 

The results obtained in the 90 to 100 km 
region in the original analysis, reference 3, are 
also not correct due to Inaccurate values of Cd. 

The recent work on CD reported in reference 6 shows 
that the values of Cq used in the original 
analysis were about lOZ higher than Chose measured 
in the wind tunnel. Therefore, somewhat higher 
values of density would be obtained in the 90 to 
100 km region but not to the degree indicated by 
the results shown in figure 2 . 

Above 100 km, an average lOZ higher denalty 
is found using the Cq values calculated from the 
transition flow analysis. This difference can 
partly be explained by the difference in the 
treatment of the gas surface interaction. In the 
current analyala, a value of K > 0 was used while 
in reference 2 the treatment of the gas surface 
interaction resulted in an additive term of the 
form 2>/ir/3Sy tms uaed where is the speed 
ratio molecules would have if they obtain the 
temperature of the sphere wall, Tw " IOO^’k. Values 
of K which would compare swre to the assumption 
made in reference 3 were attempted but the results 
for C 3 obtained for higher valuea of K were not as 
consistent as the ones obtained with K ■ 0. In 
some cases, no solution for C 3 could be found for 
R greater than zero. These results tend to 



Indicate that K •• 0 la a proper choice but aora 
conclualva evidence la needed In order to fix thla 
laportant parameter. 

Concluelona 

The aethod of analyala reported here haa 
deaKinatrated potential for application In high 
altitude falling aphere experlaenta. The valuea 
of C 3 obtained represent one of the flrat experi- 
mental measurements of this quantity under high 
altitude conditions. More accurate valuea of C 3 
are needed to remove this unVutown In the analysis 
of falling sphere experiments could then be 
designed to measure srlll other unknown quantities 
such as temperature. Dual falling sphere meaeure- 
ments of both temperature and density at high 
altitudes have been propoeed and appear to be an 
experiment well «K>rth performing. 
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Figure 3. Departure fro* atandard ataoaphere of 
of the aean and atandard dev lat Iona of 
13 denaity neaBurenenta deduced In the 
original experlaenta at Kwajaleln 
laland (plot la copied fro* reference 2). 


Figure 2. Average of ratio of denaity coaputed 
ualng equation (19) and denaity 
obtained In original experiment for 
aoundlnga 8, 13, and 14 with 
C3 - 4 X 10* a*/kg. 


49 



CHAPTER V 


A SPHERE DRAG BRIDGING RELATIONSHIP IN THE LOW 
MACH NUMBER TRANSITION REGIME 
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Introduction 

Transition flow is defined as the flow regime in which the mean 
free path of the gas molecules is of the same order as a typical dimension 
of the body under consideration.^ At the boundaries of the transition 
regime are the free molecule regime (mean free path much larger than the 
characteristic dimension) and the continuum regime (mean free path imich 
less than the characteristic dimension). The drag coefficient of a sphere 
is known to change markedly in the transition flow regime between the 
limits of free molecule and continuum flow and is a strong function of 
Mach Number. The drag coefficient in the free molec ile regime approaches 
infinity at zero Mach Number and approaches a value near 2 at hypersonic 
Mach numbers. At the continuum limit, on the other hand, the sph- re drag 
coefficient has a more complex nature which is known to depend on the 
Reynold's Number and the turbulence or lack of turbulence in the flow. 

In this work, however, the high density boundary of the transition regime 

4 

will be assumed to be at Re ~ 10 for which the drag coefficient has 

00 * 

a value of near O.A at subsonic velocities, increases in the transonic 

regime to a value of near 1.0 and approaches 0.92 at hypersonic velocities. 

Research supported under National Aeronautics and Space Administration 

Contract NAS8-28248, Marshall Space Flight Center, Huntsville, Ala. 

Assistant Research Professor, Mechanical Engineering Dept., Member AIAA. 

51 



4 

The continuum limit of the transition regime was taken as Re * 10 for 

00 

two reasons: (1) the sphere drag data employed in this work all corre- 
sponds to Re < 10^ and (2) sphere drag variations which occur above 

Re « 10^ are more clearly correlated with continuum parameters (Re and 
00 

turbulence) rather than what is normally considered transition flow 
parameters (Kn and surface to gas temperature ratio) . 

Due to the lack of adequate theory of the aerodynamics in the 
transitional regime, analytic determination of the sphere drag coefficient 
is usually made through semi empirical relations which are based on near 
free molecule flow theory and experimental results. These formulas are 
called bridging relationships, a number of which are reviewed in refer- 
ences 2, 3 and 4. The accuracy of a bridging relationship may be de- 
termined by comparison with experiment and most formula have at least 
one free parameter in order to obtain a best fit with given data. As 
discussed in reference 4, it is found, however, that available bridging 
relationships are typically accurate only over a limited range in 
Mach Number and Knudsen Number. The purpose of this paper is to report 
on a bridging formula which, with three free parameters, was found to 
predict to about 67. accuracy the sphere drag results obtained by the 
ballistic range method by Bailey and Hiatt. ^ Although Bailey and Hiatt 
provide plots of curve fits of drag coefficient as a function of Reynold's 
Number to a claimed accuracy of ± 27., the analytic results obtained here 
are provided as a function of Knudsen Number and have the advantage of 
allowing for ready interpolation as a function of Mach Number and Knudsen 
Number, Finally, it should be noted that the results obtained here are 
applicable only to fitting the Bailey and Hiatt measurements which are 

somewhat unique for sphere drag data in that the sphere surface temperature 
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was equal to the gas temperature (T^/T^ “ !)• The results then have 
application, for example, to falling sphere data analysis where T /T m 1 
but would not be applicable to wind tunnel data where typically T^/T^ » 1. 

The Bridging Relationship 

The bridging relationship used in this work is a modification of 
that developed by Matting^ and also given by Rott and Whlttenburg.^ 

Using a first collision, two fluid flow approximation. Matting obtained 
the result which can be written as 




c ^ (% 


- E/Kn 


where C is the drag coefficient, is the continuum drag coefficient, 

c 

C is the free molecule drag coefficient, Kn is the Knudsen Number 
^FM 

(defined as Kn » mean free path/sphere diameter), and E is the free 
parameter. Equation 1 is seen to provide the correct limits as Kn is 
allowed to vary. At the free molecule flow limit, Kn -♦ », which gives 

4 

C., = . At the continuum limit which in this work is taken as Re * 10 , 

D D ® 

FM 

Kn -* 0, which gives . The limits are approached asymptotically 

which is what is observed experimentally. As will be shown, however. 
Equation 1 is not found to accurately predict the Cj^ variation in the 
low Mach Number transition regime for any value of the free parameter E. 
Equation 1 is found to predict a much steeper variation in than what 
is observed in recent experimental results of Bailey and Hiatt. ^ 

In order to correct the failure of Equation 1, a second free 
parameter was introduced which was found to improve the accuracy con- 
siderably. The new form of the bridging relation is given by 




- E/(Kn) 


53 



where x is the new free parameter introduced here. By raising Kn to a 
power, the steepness of the variation of with respect to Kn may be 
controlled, thereby better fitting the experimental results. 

Method of Determination of Free Parameters 
The values of free parameters are determined from a best fit to 
experimental data. The experimental data used here is that reported by 
Bailey and Hiatt^ which are obtained by the ballistic range method for 
which T /T * 1 and covers a range in Mach Numbers from 0.1 to 6.0 and 

W 00 

a range in Reynold's Numbers from 15 to 50,000. Due to a lack of coverage 

in the transition regime at the lowest Mach Numbers, the data used in 

this work is limited to between M .72 and M fa 6.0. Also, since only 

00 00 

4 5 

6 of the 356 data points in the range have 10 < Re^ < 10 , a value of 

4 

Re « 10 has been used as the continuum limit. This range in flow 

'X) 

parameters is of particular Interest to falling sphere data analysis and 

also has application to satellite reentry and sounding rocket trajectories. 

Bailey and Hiatt provide tables of the experimental results 

arranged in groups of approximately the same Mach Number. For example, 

32 measurements of sphere drag coefficients were made for 1.45 < M^ < 1.65 

and 36 measurements made for 2.8<M <3,2 (see Table I for complete 

— 00 * 

list). For each measurement the values of M , Re , and C are given from 

GO 00 U 

which a Knudsen Number can be derived using^ 



where I is the mean free path, d is the sphere diameter, and y « 1.4. 

The continuum and free molecule drag coefficients are assumed to 

g 

be functions of Mach Number only. The free molecule drag coefficient 


is given by 


54 



*7 ■ ?) •"*•(+ • i)-=- 

(A) 

where erf is the error function and where S Is the speed ratio given 
by S ■ V/ ^2 RT *, which can also be expressed in terms of Mach Number 
as S “ M ^y /2 ’. The quantity K in Equation A is a factor of order unity 
dependent on the gas surface interaction. 

Since the evaluation of Equation A is complicated by the presence 
of the error function, a useful expansion of Equation A for low and high 
values of S was employed in the analysis. The results of the expansion 


are 


FM 


(S < 1.25) - 


[ 3 S 15 210 ^ J 


(S 


FM 




> 1.25) - 2 


(1 + K) 


(5) 


( 6 ) 


which are accurate to better than .1% with respect to Equation A. 

The continuum values of drag coefficient were obtained also from 
Bailey and Hiatt using values of C versus M for Re « 10,000 which, for 

D 00 

>A 

the Mach numbers of interest here, correspond to Kn 10 . The ex- 

pression employed is given by 

Cp (M > 1.0) * .92 + .166/M - . 366/M^ (7) 

which is found to give an accuracy of at least 57. with respect to the 
experimental results. For Mach Numbers below 1.0, the following values 
were used 

Cjj (.1 < M< 1.0) - .A03 - .1A2 M^ + .A59M^ (8) 

c 
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Using the relations given above, the data in a given Mach Number 
group were used to find the best values of x and E in a least squares 
Sinse. The least squares equation was written and the partial derivatives 
with respect to x and E were found. A computer program was developed to 
find the values of x and E which made the derivatives zero and thereby 
made the error a minimum. 

Results and Discussion 

For each Mach Number set tested, a root-mean-square (rms) value 
was computed and used as a measure of the accuracy of the fit for that 
set of data. The value of K required for the free molecule drag co- 
efficient value was found to have influence on the results obtained. 

The RMS value for a given Mach Number was found to be improved if the 
value of K was taken to be a small negative number. Therefore, K becomes 
a third free parameter of the fitting process. 

The least square results showed that x and K have a Mach Number 
dependence while E is nearly constant. The x and K values were found to 
be nearly linear in Mach Number and approximated by the following rela- 
tionships , 

X - .399 + .016 M (9) 

K - -.002438 - .01842697 M (10) 

and the average value of E was found 

E » .212 (11) 

A set of X, K, E values are thus obtained over the full range of Mach 
Numbers under consideration. Equations 9, 10, and 11 were employed in 
Equation 2 and the results compared to the ARO drag data. The RMS values 
that resulted are given in Table I along with the number of data points 
and the values of x, K, and E. The average RMS value of the 16 sets of 
data tested in Table I is .059. 
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Figure 1 shows a plot of vs Kn for six of the seventeen sets 

of data employed in the analysis. The bridging relationship is shown 

as the solid lines which are calculated based on the midrange Mach 

number of a given set of data. The figure illustrates the success of 

the bridging relationship in the transition regime and shows that much 

of the 6% rms error can be traced to scatter in the data which is nearly 

± 107„ at some values of Kn. One failure worth noting, however, is the 

tendency of the bridging relationship to underestimate the Cjj value in 
-3 -2 

the 10 < Kn < 10 range by about 5% in some cases. This is likely 

a slip flow influence which has not been taken into consideration in 
this work. 

Table I 

RMS Values of Curve Fit to ARO Data Using Equations 9, 

10, and 11 in Equation 2 


Midrange Mach Number of 


Mach Number 

Number Range Data 

Points 

X 

E 

K 

RMS 

.72 

.70 

- 

.74 

7 

.4105 

.212 

-.016 

.013 

.81 

.79 

- 

.83 

9 

.4120 

.212 

-.017 

.076 

.915 

00 

00 

• 

- 

.95 

20 

.4136 

.212 

-.019 

.065 

.965 

.96 

- 

.97 

10 

.4144 

.212 

-.020 

.123 

.989 

.98 

- 

.998 

8 

.4148 

.212 

-.021 

.082 

1.135 

1.08 

- 

1.19 

24 

.4172 

.212 

-.023 

.036 

1.25 

1.2 

- 

1.3 

15 

.4190 

.212 

-.025 

.055 

1.375 

1.3 

- 

1.45 

28 

.4210 

.212 

-.028 

.075 

1.55 

1.A5 

- 

1.65 

32 

.4238 

.212 

-.031 

.050 

1.75 

1.65 

- 

1.85 

23 

.4270 

.212 

-.035 

.051 

2.05 

1.9 

- 

2.2 

30 

.4318 

.212 

-.040 

.060 

2.55 

2.4 

- 

CM 

18 

.4398 

.212 

-.049 

.059 

3.00 

2.8 

- 

3.2 

36 

.4470 

.212 

-.058 

.054 

A. 00 

3.8 

- 

4.2 

28 

.4630 

.212 

-.076 

.060 

5.00 

4.8 

- 

5.2 

40 

.4790 

.212 

-.095 

.041 

6.00 

5.8 

- 

6.2 

33 

.4950 

.212 

-.113 

.042 
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Conclusions 


A sphere drag bridging relationship has been developed for the 
low Mach Number transition flow regime which fits recent experimental 
results to an accuracy of about 6^ rms. The experimental results used 
were exclusively those reported by Bailey eud H’att in tests ran at ARO 
and reported in March 1971 in which the sphere surface temperature was 
equal to the gas temperature. The results of this work should have 
application, for example, to the analysis of falling sphere data in 
which T /T w 1 . 

V 00 

Due to the unique nature of the Bailey and Hiatt data (i.e., 

T /T = 1), it is of interest to examine the conclusions these data 

W 00 

indicate concerning the nature of the transition flow regime. Using 
the parameters found in fitting these data, Eq. 2 was plotted, vs M, 
(Fig. 2) for 0 < Kn < 00 and .1 < M < 6 for constant values of Kn. Since 
the highest Kn tested by Bailey and Hiatt was 10 ^ and since there was 
little transition flow data below M “ 1, the curves for Kn > 10 ^ and all 
the curves for values of Kn for M < 1 are extrapolations. The results 
obtained, however, point to two important conclusions: (1) the width of 
the transition regime in terms of Knudsen Number is wider than usually 
assumed and (2) the free molecule drag coefficient implied by the results 
is less than usually assumed. 

Since the value of x (which is a measure of the slope of the Cj^ 

vs. Kn curve in the transition regime) was found to be about .45, the 

width of the transition regime is increased over that obtained using the 

Matting relation which has x ■ 1 (Equation 1). This conclusion is 

illustrated by substituting into Equation 2 the value Kn = 5 which is 

the usually assumed upper limit of the transition regime.^ The results 

58 



at M ■ 6, Kn * 5 give C /C ■ .956 which shows that the free molecule 

° *^fm 

limit has not been reached at this value of Knudsen Number. A value 

of Cjj/Cp ■ .99 is found to be reached for M ■■ 6 at a value of Kn ■ 
fm 

100 . 

The free molecule limit was found to be of importance in the de- 
velopment of the bridging relationship since K had to be adjusted to neg- 
ative values in order to obtain an accurate fit. This implies a free 

molecule drag coefficient less than tt/o for the high values of Mach 

9 

Number tested. Results from other experimenters show that the drag 
coefficient in the free molecule limit is greater than 2 at Mach Numbers 
from A to 6. This departure from the results of others is likely ex- 
plained in that much of the sphere drag data at high Knudsen Numbers and 
high Mach Numbers used by others are obtained in low density windtunnels 
whereas the experimental data employed in this study was obtained ex- 
clusively from ballistic range data. The higher surface- to-gas-tempera- 
ture ratios that occur in windtunnels cause higher free molecule drag 

due to the energetic reflection of molecules at the surface. This effect 

9 10 

is clearly shown in recent sphere drag experiments ’ in which C 

°fm 

is found to decrease as T /T is decreased. In fact, an extrapolation 

W 00 

of the data of reference 10 down to T /T » 1 indicates a C of 2 or less. 

W 00 D 

A free molecule sphere drag coefficient of less than 2 requires that the 
gas surface interaction be non uniform on the surface of the sphere. 

This possibility is discussed by Cook^^ in which he proposes that the 
lowest possible free molecule drag of a sphere is 1.5. The lowest value 
implied in this work is 1.84A at M = 6 where K “ -.113. 

The above conclusions are tentative but point to a need for free 
molecule experimental data in the low Mach Number regime at surface 
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temperatures close to ambient. Although the Bailey and Hiatt data do 
not reach free molecule conditions, the extrapolation discussed above 
show that the data lead to different conclusions than obtained from 
T /T ^1 experiments, 

W CO 
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FIGURE CAPTIONS 


Figure 1. Comparison with Equation 2 of Sphere Drag Coefficient Data 
(Ref. 5) as a Function of Free Stream Knudsen Number for 
Constant Values of Mach Number. 

Figure 2. Sphere Drag Coefficient as a Function of Free Stream Mach 
Number for Constant Values of Knudsen Number. 
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CHAPTER VI 


IMPROVEMENTS IN FALLING SPHERE DATA ANALYSIS IN THE 80 TO 120 KM REGION 


64 



Improvements in Falling Sphere Data Analysis in the 80 to 120 Km Region 

by 

Gerald R. Karr and Robert E. Smith 

The analysis of falling sphere drag data is a principle means of 
density and temperature determination in the 80 to 120 Km region of the 
earth's atmosphere. This important method of atmospheric probing was 
reported by Bartman, et al^ in 1956 and has found wide use in upper atmo- 
spheric research. The method is particularly useful in the 80 to 120 Km 
region which is above the altitude capability of most aircraft and balloon 
probes but below normal satellite altitudes. 

In a typical falling sphere experiment, a sphere is ejected from a 
sounding rocket and the trajectory of the sphere is measured as it falls 
through the atmosphere. The trajectory information is analyzed to de- 
termine the velocity and acceleration as a function nf altitude. This 
information is then used in the drag equation 

Drag - p Cjj A (1) 

where density, p , is obtained for a given value of drag coefficient, C^, 
and sphere crosssection area A. 

Temperature values are obtained from the density measurement using 
the hydrostatic equation 

dp » - g p d z (2) 

where g is the acceleration of gravity, z is the altitude, and p is the 
pressure. Equation 2 may be integrated between any two limits in altitude 
to give the pressure at z 
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(3) 


z 

p(z) «- y* gpdz+p (z^) 
z 

o 

The temperature at z is obtained by substituting the values of p and p 
into the equation of state, 

p(z) - p (z) (R/W(z)) T(z) (4) 

where R is the universal gas constant, W is the molecular weight and T 
is the temperature. 

The temperature determination is seen to require knowledge of the 
pressure or equivalent temperature at some reference altitude (p(z^) in 
Equation 3). In practice, the integration of Equation 3 is taken in the 
negative direction, from the point of highest ascent down to the lower 
altitudes. The temperature at z^ is usually obtained from some atmospheric 
model. The error caused by possible incorrect temperature selection is 
minimized by the practice of downward integration since the term 
in Equation 3 becomes small in comparison to the first term as the inte- 
gration proceeds to the lower altitudes. Thus, the effect of error in the 
initial temperature selection should become unimportant at one to two scale 
heights below the initial selection point. This was verified in the present 
study by selecting temperatures from 0°K to 1000°K at 120 Km. The effect 
on results below 100 Km due to such wide choices in temperature at 120 Km 
was found to be negligible ( < 1%). 

Another source of error in the falling sphere method is the value of 
Cp used in the drag equation. Equation 1. The typical trajectory of falling 
spheres in the 80 to 120 Km region is found to correspond to an aerodynamic 

flight regime for which sphere drag coefficients have been highly uncertain. 

-5 -8 3 

Due to the low densities (10 to 10 Kg/m ) and the low Mach numbers 
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(1 to 5) experienced by falling spheres in this region, the aerodynamic 
flight regime Is classified as the transition flow regime (Ref. 2). The 
mean free paths of gas molecules In this region of the atmosphere vary 

_3 

between 10 m to lOro which means that the flow Is too dense to be considered 
free molecule but two rarefied to be considered continuum. While the free 
molecule and continuum values of sphere drag coefficients have long been 
well known, the values of in the transition flow region have only re- 
cently been measured to adequate accuracy. 

The purpose of this paper Is to report on Improvements made In the 
falling sphere method of analysis. The most Important of which is an Im- 
proved relationship for the sphere drag coefficient which has estimated 
accuracy of at least 5% over the range applicable to falling sphere tra- 
jectories. The second Improvement in the analysis is to employ both the 
ascent and descent trajectory data In the data reduction. These improve- 
ments are applied to published fa. ling sphere data and comparison Is made 
with the results of previous analysis. 

Drag Coefficient Relationship 

3 

Recent sphere drag experiments reported by Bailey and Hiatt have 

provided and Improved the drag coefficient Information in the transition 

regime. The experiments were made in a Mach number and Reynold's number 

range of particular application to falling sphere data analysis (0.1 < 

M < 6.2 and 20 < R < 10^). These data were obtained in ballistic range 
00 e 

in which the temperature of the gas and the sphere temperatures were 
approximately the same. The data are obtained in a consistent manner by 
one group of experimenters over the complete range of flow parameters 
applicable to the falling sphere flight regime. In view of the Importance 
of this single source of drag data, a curve fitting relationship employing 
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4 

this data exclusively has been prepared by Karr. The sphere drag bridging 
relationship thus developed Is given by 


C 


D 



+ e 


-E/a (Kn)^ 




(5) 


g 

where a = .499 and where C., and C are the continuum and free 

" “c ®FM 

molecule drag coefficient, respectively. The quantities E and x are 

parameters of the curve fit and Kn Is the free stream Knudsen number. 

The free molecule and continuum drag coefficients may be written as 

4 

functions of speed ratio, S , and Mach number, M, given by 

.3 


C_ (S < 1.25) - 

FM 3 n 




( 6 ) 


C (S > 1.25) = 2 

FM 


{ 


1 + 




(1 + K) 


(7) 


Cjj ( M > 1.0) » .92 + .166/M - 366/M' 


( 8 ) 


Cp (.1 < M < 1.0) » .403 - .142M^ + .459M^ 
c 


(9) 


where 

S » V/^2 RT ’ and M = V/ Ej'’ (10) 

^ W W 

where y Is the ratio of specif 1.: heats of the gas, and (1 + K) In 
Equations 6 and 7 Is a factor or order unity which Is related to the gas- 
surface Interaction. The least squares curve fit of the above relations 
to the data of Bailey and Hiatt has revealed that Equation 5 will fit the 
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data to an accuracy of about 57* (based on RMS value) for the following 
values of the parameters 


X * 

.399 + .016 M 

(’1) 

K “ 

-.002438 - .01842697 M 

(12) 

E * 

.212 

(13) 


The use of Equations 6 through 13 In the sphere drag bridging rela- 
tion given In Equation 5 provides values of drag coefficients which are 
based on the most accurate and applicable Information now available. In 
order to use this relationship In the drag equation, values must be 
available for the velocity V, temperature T, molecular weight W, and the 
Knudsen number, Kn. The velocity Is a measured quantity while values for 
the others are obtained as discussed In the following. 

Knudsen Number 

Knudsen number Is defined as the ratio of near free path, to the 
characteristic length of the object. We take the characteristic length 
to be the sphere diameter, d . Therefore, 


Kn - X / d (14) 

The mean free path Is Inversly proportional to number density, n, for a 
simple gas (c.f. Ref. 5) and given by 



where o Is the diameter of the gas molecules. Using the average molecular 
weight, we write the mean free path as a function of density 


I 


(a n 




(16) 
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where N. is the Avogadro constant. Using a value of a * 3.72 x 10 m 
A 

which is representative of air^, we obtain 

Kn » 2.7 X 10’^ (17) 

pd 

3 

for 0 in units of Kg/m and d in meters. 

In developing Equation 17, we have neglected the small effect tempera- 
ture has on the value of a and we have assumed the gas molecules each 
have a mass corresponding to the average molecular weight of the gas. 

Molecular Weight 

A value of molecular weight is required for the Mach number, speed 
ratio, and Knudsen number expressions. (We assume y “ independent 

of the molecular weight). The mean molecular weight is known to be a con- 
stant of 28.964 up to about 90 Km altitude. At that altitude, the heavier 
molecules begin to settle out causing a drop in molecular weight. Models 
of this molecular weight with altitude are provided in Ref. 6. Representa- 
tive values of the variation is given by the following equation 

W = 24.68 + .1235 z - .000874 (18) 

where z is the altitude. Equation 18 gives values of molecular weight 
with an accuracy of better than 1% in comparison with a nominal spring/ 
fall values given in Ref. 6. 

Temperature and Density Iteration 

The temperature at a given altitude is determined as stated in the 
introduction, by downward integration of the density profile. However, 
since the drag coefficient is a function of temperature (i.e., Mach number 
and speed ratio are inversely proportional to the square root of temperature). 
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values of temperature must be available before determination of density 
is made using the drag equation. Thus, we are lead to an iterative 
procedure for finding both density and temperature. The procedure is 
begun by assuming some initial temperature profile from which, given the 
unusual velocities, the Mach number and speed ratio are found. Once a 
density profile is obtained (to be discussed in the next paragraph) using 
the assumed temperature profile, a new temperature profile may be con- 
structed. This process is repeated until the values of both temperature 
and density no longer change beyond a specified error limit. 

Since the drag coefficient has been written as a function of density, 
the drag equation becomes of the following form 

Drag = p A C + (C - C ) (19) 

c FM c 

where b is defined through Equation 17 and 5. 

b - E/(2.7 X 10'^ W/d)^ (20) 

Equation 19 is a nonlinear equation in the unknown density and must be 
solved using numerical procedures. 

Solving for density from the drag equation of the form given in Equation 
19 has two advantages. First, as already pointed out, the relationship 
for Cjj represents recent, accurate drag measurements. Second, by solving 
Equation 19 for density directly, we eliminate the uncertainty of choosing 
a drag coefficient from a set of values tabulated as a function of Mach 
number and Reynold's number. Falling sphere experimenters effectively 
solve an equation like Equation 19 by iterating between tabulated values 
and using the simple drag equation. The solution of Equation 19 is obtained 
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much faster and the solution is likely more accurate than that obtained 
in such iterative procedures. 

Application of the Proposed Method 

The procedure outlined above was applied to sets of falling sphere 
data reported by Peterson, et al^> for measurements made over Kwajalein 
during 1963 and 1964. One of the reasons for using this data to illustrate 
the application of the proposed method of analysis is the importance of 
the results which were obtained in the original experiment. This set of 
data is of particular interest since it forms the basis for the density 

and temperature model in the 80 to 120 Km region for what is labeled 

o 6 

"15 N, Annual" in the U. S. Standard Atmosphere Supplements, 1966. A 

summary of the density results obtained by Peterson, et al? are given in 
Figure 1 which is a copy of Figure 2.20 of Reference 6. 

A second reason for using this particular set of data was the avail- 
ability of the information required in our analysis. Through one of the 

g 

original experimenters, K. McWatters, we obtained the detailed computer 

output (a sample of which for one flight is given in Reference 7) for the 

falling sphere measurements made by the group. The data of interest is 

the density, p^, temperature T^, drag coefficient Reynold's number 

Mach number M , and velocity V , as a function of altitude, z , which 
on n 

resulted from their analysis of the falling sphere trajectory. Of these 
quantities, only velocity and altitude were measured where the subscript 
n is used to indicate this. The remaining quantities were deduced and we 
use the subscript o to indicate this. The measured drag force, however, 
can be found by taking the deduced density and the drag coefficient em- 
ployed in obtaining that density, and substitute these values into the drag 
equation. Therefore, 
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(Drag)^ 



A 

n 


where the area Is based on 0.66 m sphere diameter. Thus, in this manner, 
the values of altitude, velocity, and drag force are obtained as measured 
in the original experiments. The accuracy of these measurements is ex- 
pected to be at least 57,? 

In order to begin the iteration for density and temperature, the data 
analysis requires that an initial temperature profile be given. Of all the 
temperatures given in this initial profile, the only one of importance is 
the one given at the altitude at which the downward integration is begun. 
For the Kwajalein data, the highest altitude for which data are available 
is 120 Km since data above 120 was considered too inaccurate. In our 
analysis we have taken the temperature at 120 Km to be 460'\ which is based 
on the 15°N Annual model given in Reference 6. For comparison, Peterson, 
et al? used a temperature value of 36l°K at 120 Km which is based on the 

9 

1962 U. S. Standard Atmosphere. The higher temperature was chosen in this 
work in view of the fact that the results obtained in the original analysis 
revealed that the temperatures were generally higher than the USSA 1962 at 
the high altitudes. The effect the choice of a 100°K higher temperature has 
on the final results will be discussed more in the conclusions. 


Temperature Interpolation for the Descent Trajectory 

At this point in the discussion, we describe another improvement to 
the data analysis which is applicable to the Kwajalein falling sphere data. 
Much of the data obtained over Kwajalein consists of two trajectories; 
ascent and descent. Soon after ejection of the sphere, the radar tracks 
the sphere during its ascent to apogee. This tracking begins usually near 
100 Km and data is then terminated at 120 Km for accuracy reasons as dis- 
cussed above. During the passage through apogee, the radars "lose" the 
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sphere. Also because of smoothness requirements, data processing during 

(See Fig. 2) 

descent generally does not start again until about 100 Km. In the orig- 
inal analysis of this data by Peterson, et al7, the ascent and descent 
trajectories were analyzed separately. That is, the downward integration 
required for temperature determination was begun again at the top of the 
descent trajectory with a temperature choice based on USSA 1962. The 
temperature thus chosen at the top of the descent trajectory is potentially 
erroneous, the effect of which will propagate two scale height down into 
the descent trajectory. We have developed a more accurate procedure for 
choosing the temperature value at the top of the descent trajectory which 
employs the ascent temperature values and an isentropic relationship. The 
ascent temperature values in the 100 Km region are considered to be mor 
accurate than the higher altitude values since this corresponds to two 
scale heights below the 120 Km altitude where the temperature is taken 
arbitrarily to be A60°K. 

The temperature determination procedure employed in the present 
analysis is as follows: 

Step 1. Density profiles are determined for both the ascent and de- 
scent trajectories, p. (z) and p_.(z), respectively. The 15° N Annual model 

u 

is used to give an initial temperature profile. 

Step 2. A new ascent temperature profile, determined for 

the ascent data using the downward integration method. 

Step 3. For the region of altitude for which both ascent and descent 
data points are available, an isentropic relationship is employed to find 
the temperature for the descent trajectory points, Tp(z), given by 

Y - 1 

T„(z) - T^(z) (opCO / p^(0) 
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where Y the ratio of specific heats. In using this relationship, we 
are assuming that the densities which have been measured at the same 
altitude by the drag method are different due to an Isentroplc process. 

The sphere trajectory is such that at about 100 Km altitude, the ascent 
and descent trajectories are about 80 Km apart. Wave motion in this 
region of the atmosphere could then account for the differences In density 
which are seen. Over the length and time scales of Interest here, the 
assumption that such processes are Isentroplc Is reasonable. 

Step 4. For the remaining altitude points In the descent tra- 
jectory, the temperatures are determined using the downward Integration 
method. 

Step 5. Based on the new temperature values, calculate new 
values of Mach number and drag coefficient to be employed In the 
determination of new density profiles. That Is, start again at Step 
1 above and continue the iteration procedure until convergence is 
reached in both temperature and density results. 

The above procedure requires at Step 3 that at least one data 
point of the descent trajectory be at the same, or nearly the same, 
altitude as a point on the ascent trajectory. This requirement Is met 
for six of the thirteen flights made over Kwajalein and reported in 
reference 7. Table I gives some of the characteristics of the flights 
studied In this work including the range of overlap for the six flights. 
Two of the six flights had only one point in common while flight # 13 
was found to have 12 data points In common (dat.. Is provided approximately 
every kilometer, on the kilometer). 


75 



Results 


Density and temperature profiles were obtained by the methods de- 
scribed above for the six falling sphere flights over Kwajaleln which 
had overlap In altitude coverage of the ascent and descent trajectories. 

The density profiles were compared with that published In the U. S. 

9 

Standard Atmosphere 1962. The ratio of the density found In this work 
to the USSA 1962 densities Is shown for each of the six flights In Figures 
3 through 8. Also contained In these figures are the temperature pro- 
files obtained in the analysis. Table 1 lists the sounding number which 
was designated by the experimenters, the time and date of the flight, 
and the altitude range covered on the ascent and descent trajectories for 
the six flights analyzed. 

The mean and standard deviation of the density and temperature pro- 
file ratios (with respect to the models Indicated) are plotted In Figures 
9 and lO, respectively. For the altitudes outside the region of over- 
lap, six values were used In construction of the mean and standard de- 
viation. For some altitudes within the overlap region, as many as 12 values 
were used. 
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Discussion of Results 


The summary of density results given In Figure 9 show a signifi- 
cant departure from that obtained In the original analysis of Peterson, 
et al. 1966. This departure is found to be primarily because of the In 
complete Information on drag coefficients available to the original 
experimenter at that time. The drag coefficient values employed in 
the present analysis are considerably more accurate and reveal a point 
of maximum departure from the 1962 standard at 105 Km rather than at 
92 Km. 

The temperature results given In Figure 10 show better agreement 
with the original analysis than does density. The departure at 120 Km 
Is artificial since the choice of a temperature value at 120 Km Is 
completely arbitrary. The departure In temperature between the present 
and original analysis at high altitudes Is not the cause for the de- 
parture In density at these altitudes. In fact. If the value of 
temperature at 120 Km used by Peterson, et al. were to be employed In 
the present analysis, the departure In density results would be even 
greater than given in Figure 9. The colder temperatures used by 
Peterson, et al. would result In lower drag coefficients and therefore 
even higher densities would be obtained from the drag equation. The 
results we have obtained seem to clearly point towards the higher 
values of temperature at altitudes near 120 Km. 
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Table 1 


Pitura IhMibar 

Soundint Number 

Time 

(GKT) 

& Date of Flight 

Range of 
Ascent 
Altitude 

Range of 
Descent 
Altitude 

Overlap 

Range 

3 

2 

0257 

March 29, 1963 

100-120 

102-80 

100-102 

4 

3 

0328 

June 18, 1963 

99-120 

102-80 

99-102 

5 

8 

1458 

Nov. 14, 1963 

104-120 

109-80 

104-109 

6 

12 

1820 

March 13, 1964 

99-120 

104-80 

99-104 

7 

13 

1125 

May 12, 1964 

96-120 

107-80 

96-107 

8 

14 

0101 

June 17, 1964 

102-120 

102-80 

120 
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Departure (percent) 


Fig. 1. Departure from standard atmosphere of the mean and standard 

deviations of 13 density measurements deduced in the original 
experiments at Kwajalein Island (plot is copied from Ref. 2). 



Fig. 2. Velocity-altitude history for falling sphere sounding 13 
in which ascent and descent trajectory overlap. 
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Temperature , K Density Ratio, p/p (ussa.issz) 


Figure 3. Temperature and Density Profiles Obtained fron 
Sounding if 2. 



Temperature, K Density Ratio, /='// 5 {ussa.i 962 ) 


Figure 4. Temperature and Density Profiles Obtained from 
Sounding 0 3. 

80 



Temperature, K Density Ratio./y/>(ussA.(9€2) 


Figure 5. Temperature and Density Profiles Obtained from 
Sounding it 8. 



Temperature , K Density Ratio, />//> (ussA.isez) 


Figure 6. Temperature and Density Profiles Obtained from 
Sounding if 12. 
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Figure 7. Temperature and Density Profiles Obtained from 
Sounding # 13. 
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Figure 8. Temperature and Density Profiles Obtained from 
Sounding // lA. 


82 



ALTITUDE, Km 



DENSITY RATIO, y’/y’CUSSA. 1962) 


Figure 9. Mean and standard deviation of density results obtained 
In the present analysis for six falling sphere flights 
over Kwajaleln. Mean of results from original analysis 
given for comparison. 
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TEMPERATURE/tEMPERATURE(USSA.I962) 


Figure 10. Mean and standard deviation of temperature results 
obtained in the piesent analysis for six falling 
sphere flights over Kwajalein. Mean of results from 
original analysis given for comparison. 
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FREE MOLECULE DRAG AT SPEED RATIO LESS THAN UNITY 
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Free Molecule Drag at Speed Ratio Less Than Unity 

by 

Gerald R. Karr 

Assistant Research Professor 
Mechanical Engineering Department 
The University of Alabama in Huntsville 

This work concerns results obtained in the calculation of the 
drag force acting on objects which (1) have dimensions much less than 
the mean free path of the gas (Knudsen number much less than one) and 
(2) have velocity with respect to the gas which is much less than the 
thermal velocity of the gas (Speed ratio, S, much less than one). These 
conditions are characteristic of those experienced by certain aerosols 
and Brownian type particles which are suspended in a gas or which are 
forced to travel through a gas such as in separation processes. 

The purpose of this work is to investigate the influence of two 

factors which enter the calculation of the free molecule drag force; 

(1) the shape of the object and (2) the gas surface interaction. Past 
1 2 3 

investigations ’ * have concentrated on the second of these factors 
while assuming the particles are perfectly spherical in shape. The 
gas surface interaction has generally been taken to be composed of a 
specular fraction and a second fraction which is purely diffusive (see, 
for example. Ref. 4) or a modified diffusive such as the elastic- 
diffusive reflection employed in Ref. 3. The form of the diffusive 
fraction of the reflected molecules has received considerable attention 
because comparison of sphere drag calculations with experiment lead to 
the conclusion that all the molecules must be diffusively reflected in 

order to explain the observed high drag coefficients at low speed ratios. 
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In the work reported in this paper, the drag coefficients for both 
specular and diffuse reflection were obtained for non»spherical objects 
in order to investigate the influence of shape at very low speed ratios. 

The drag coefficients were obtained by employing the expression 
for force acting on an element of surface in a free molecule flow valid 
at any speed ratio and taking the component of that force which acts 
in the direction of the velocity vector. The drag component of force 
was divided into a part due to the momentum of impinging molecules, 

D^, and a part produced by the reaction force of the molecules leaving 
the surface, D^. Expressions for both specular D^(specular) and 
diffusive D^(diffusive) reflections were developed. The total drag 
force is then given by D = D. + D for an element of surface at any 
angle with respect to the flow and any speed ratio. The expressions 
were then integrated over the surfaces of various shapes including 
oblate spheroids, cylinders with flat and spherical ends, and cones 
with flat ends. The drag force coefficients were obtained for these 
objects at low speed ratios for both specular and diffusive reflection. 
The total drag force acting on an object is written 

D « 0, (1 + ) 

where the ratio represents the effect of the reflection normalized 

with respect to the incident contribution. For a perfect sphere, for 
example, and for S « 1 

D. (sphere) = p tt r^ ^2tt RT' V 

The attached figure presents examples of values of 0^/D^ for oblate 

spheroid shapes as a function of the minimum-to-maximum radius ratio. 
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The results show the sensitivity of the to both the reflection 

characteristic and the object shape. 

The results obtained in this work show that the magnitude of 
the force coefficient ii strongly dependent upon the shape and o»^ien- 
tation of the object for both specular and diffusive reflection. Since 
non-spherical aerosol or Brownian particles would present randon orien- 
tation with respect to the velocity vector, the force coefficient for 
a given object would be an average value. The results obtained in this 
work provide the information needed to obtain the average force co- 
efficient for various non-spherical shapes. One conclusion reached in 
this work is that both specular and diffuse reflections can produce 
high drag force coefficients at low speed ratios for non-spherical 
objects and neither should be excluded from consideration in such cases. 
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D^/D^ ratio at low speed ratio for oblate spheroid shapes as 
a function of minimum-to-maximum radius for both specular and 
and diffuse reflection niodels. 
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A 


Proposal to Develop Zero-g Brownian Motion Experiments 

by 

G. R. Karr 


Introduction 

Robert Brown In 1828 Is credited with establishing as an impor- 
tant phenomenon the observed Irregular and perpetual motion of macro- 
scopic particles suspended In gases or liquids. The theory of this 
aratlon, %4iich has become known as Brownian motion, has received the 
attention of Einstein, van Bmoluchowskl , Langevln , Uhlenbeck and 
Ornstein, Chandrasekhar, and many others. The motion received much 
early Interest because It has been established through the theory that 
Broimian motion Is direct observation evidence of the molecular state 
mf Blatter. For example, from observation of Brownian motion, one can 
-BKasure Avrogodros number and this method was In fact considered to 
provide the most accurate oieasurement of this quantity In the early 
1900* s. The theory of Brownian motion establishes that the motion Is 
described am a random process which Is found to be a major step In the 
developsent of the field of study now called stochastic processes. 

Experimental Investigations of Brownian motion has not received 
the Interest of physicists recently and modern Interest In Brownian 
motion Is primarily In the theory and Its application. However, with 
the unique environment which %d.ll be provided by the Space Shuttle, the 
possibility now exists for performing experiments Involving Brownian 
motion that could not be done In a one-g environment. It Is proposed 
here to study possible experlaients Chat may be performed In a zero-g 
environment which Involve Che observation of Brownian motion. It Is pro 

posed to assess the feasibility and Importance of such experiments. 
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The Theory of Bro%mlan Motion 


The currently accepted theory of Brownian motion was presented 
by Ornstein and Uhlenbeck in 1930. Beginning with the equation of 
motion of a B*’ownlan particle (called the Langevin equation) we have 
.2 

m — X “ -mg V + F(t) (1) 

dt 

where 3 is a drag coefficient and where F(t) is a random forcing term 
characteristic of the Brownian motion. Ornstein and Uhlenbeck obtain 
the solution for the mean square displacement of the particle given by 


g2 ^ . 2 . k _T (3t - 1 + e ■ (2) 

m 3 

which has the limits for t large compared to 3 ^ 



2 k T 
m3 


(3) 


which is the result obtained by Einstein and for t small compared to 
3 ^ the result is 



(A) 


where u^ is the initial velocity of the particle. 

Other quantities may also be calculated. For example, the mean 
square velocity of Brownian particles all starting at velocity u^ is 


2 

u 


u 

o 



g- 2 3 t 


(5) 


and the velocity distribution function of the particles is given by 
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( 6 ) 


which shows that the particles eventually reach a Maxwell-Bolttmann 
distribution. 

The above results are all for what is called a free particle. 

That is, the only forces acting on the Brownian particle are the drag 
term m0v and the random forcing term F(t). For the case of Brownian 
motion in an external force field, such as gravity, theory Is not so 
clear nor as complete as for the free Brownian motion. Uhlenbeck and 
Ornstein, for example, consider the Brownian motion of a particle which 
is bound in a harmonic force field at frequency Three cases result 

and the solutions for the mean square displacement under these conditions 
are provided. 




While Uhlenbeck and Ornstein obtained solutions for the harmonic 
forcing cases, the case for constant forcing has not been solved as 


completely. 
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For example, Wang and Uhlenbeck state In their 1945 paper that they 
have been unable to find the general solution for the constant force 
case. Solutions do exist, however, for the special case when the 
friction force Is high and the observation time Is large (l.e., for 
t » 0 ^). Chandrasekhar, for example, obtains the time varying dis- 
tribution of Brownian particles In a gravity field and Is able to show 
that the particles arrange themselves In a barometric distribution with 
respect to the gravity vector. 

Solutions for other observables of Brownian motion In a constant 
force field are apparently not available. For example, the mean square 
displacement or the velocity distribution function in a gravity field 
were not found In tuw references cited. Part of this study will be di- 
rected toward a thorough investigation of the recent literature on this 
subject. From the literature that has been searched, however, it is 
apparent that Brownian motion in a gravity field will be much more com- 
plex than the free motion and that sensitive experiments in one-g would 
be strongly Influenced by gravity. 

Zero-g Experiments 

The zero-g environment offers at least three advantages in the 
performance of Brownian motion studies. 

1. The theory of Brownian motion in zero-g is well established 
while the motion under gravity is difficult at best to interpret. 

2. Convection currents can be minimized or eliminated in zero-g, 

3. In zero-g, particles of larger sizes and masses can be em- 
ployed in Brownian motion studies. 

The purpose of the proposed work is to study possible Brownian 

motion experiments and to assess the value and feasibility of such 

experiments. While it is expected that other experiments will be proposed 
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and studied during this work, the following set of experiments have 
already received preliminary consideration and will be described 
briefly, 

A. Particulate drag In the transition regime 

As Is well known. Brownian motion is observed In both liquids 
and gases. The drag or friction coefficient, these media Is, 

however, considerably different. The friction factor in a liquid for 
example is In the fluid dynamic flow regln.e called "Stokes flov/' and 
Is described as a highly- viscous flow. For a Brownian motion In a 
gas, however, the flow Is free molecule if the mean free path of the 
molecules is large compared with the size of the particle. It is 
therefore proposed that a possible zero-g experiment is to var/ the 
properties of the suspension medium or the particles, so as to obtain 
friction factor information in the transition flow regime between the 
Stokes regime and the free molecule regime. The proposed work would 
consist of determining the range of particle and medium properties 
that would be required to probe the transition regime. Also to be studied 
is the limits Imposed on such an experiment by the one-g environment. 

The transition flow regime ha \ proved very difficult to probe in ground- 
based experiments and becomes progressively more difficult as the speed 
is reduced. The experiments proposed here would provide data at the very 
low end of the velocity spectrum, a region for which little or no data 
now exists. This possibility of increasing the range of understanding 
of fluid dynamic drag appears to be of great value. 

B. Gas surface Interaction 

Brownian motion observations provide an excellent basis for 
studying the interaction of molecules with surfaces. The motion is a 


direct consequence of the bombardment of the particle surface with 
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molecules of the surrounding medium. The free molecule regime is best 
for such observations since the motion is isolated from the effects of 
viscosity that occur in the Stokes flow regime. The proposed experi- 
ments would consist of (1) the preparation of particles of known com- 
position and surface properties, (2) the preparations of gases of known 
composition and temperature, and (3) collection of data on the resultant 
mean displacement and velocity distribution to determine the friction 
factor 0. The value of 3 determined in such an experiment can be di- 
rectly related to the effect of the gas surface Interaction. The data 
can then be correlated with respect to the gas and surface properties. 

The physics of the gas surface interaction is not well understood 
at present and experimental data of the type proposed here would be valu- 
able in identifying the important characteristics of the interaction. 

A factor of two variations in the friction coefficient, 3, is theoreti- 
cally possible due to the effect of the gas surface interaction. Ground 
based experiments such as the oil drop experiment do not have the sensitiv- 
ity needed to determine these effects. 

C. Verification of Brownian Motion Theory 

The Brownian motion theory proposed by Uhlenbeck and Ornstein is 
the currently accepted theory of the motion. The theory of Einstein and 
Smoluchowskl is found to be the limiting case of the Uhlenbeck and 
Ornstelii theory for large times (t » 3 for the free particle case. 
Edward Nelson further proposes that the Einstein-Smoluchowski theory is 
also limited to the case of large friction (3 large) and points to ex- 
perimental results of Kappler and of Barnes and Silverman to show that, 
for the harmonic forcing case, the Einstein-Smoluchowski approximation 
is invalid for the under damped case even for t » 3 For the same 

reasons as mentioned above, the presence of the one-g field requires 
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that particles be small and light. For such cases then, it is difficult 
to probe the Brownian motion for short times (t < p ^), The motion for 
short times can be used to verify the Uhlenbeck-Ornstein theory. The 
results for t » p ^ are the same for both the Uhlenbeck-Ornstein and 
the Einstein-Smoluchowski theories and thus such results can not be used 
to distinguish the theories. Due to the difficulty in making measure- 
ments in the short time on earth, there appears to be little, if any, 
experimental verification of the Uhlenbeck-Ornstein theory. It is, 
therefore, proposed that since the zero-g environment offers the 
opportunity to adjust the size of p over a wide range, a properly de- 
signed Brownian motion experiment in space would allow for the verifi- 
cation of theory of Brownian motion. The proposed work will seek to 
establish the conditions required to perform such an experiment. 

Proposed Work 

The three experiments proposed above are clearly of great value 
and preliminary study indicates that they are also feasible. During the 
proposed study, these experiments and others will be considered in detail 
to determine the value of the experiment, the reasons that zero-g are 
required, and the feasibility of performing the experiment. Extensive 
literature searches and personal contacts will be made to ascertain the 
state-of-the-art knowledge of Brownian motion theory, motion experiments, 
friction coefficients, gas-surface interaction, etc. Experimental tech- 
niques will be surveyed and recommendations made for space experiments. 

It is also expected that preliminary experimental development will be 
undertaken to test observation methods, data collection methods, and 
data analysis techniques. 
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B 

An Experiment Using the Molecular Beam Apparatus Proposed for Space Shuttle 
Title: Aerodynrmic Force Measurements in Space 

Principal Investigator: Dr, Gerald R. Karr, The University of Alabama in 

Huntsville, Huntsville, Alabama. 

Summary of Proposed Work 

The feasibility is to be determined of employing the satellite orbit 
environment in the measurement of forces resulting from the interaction of 
atmospheric gas with solid surfaces at satellite velocities. To be con- 
sidered is an experiment designed to measure the aerodynamic forces acting 
on surfaces exposed to the high-velocity, low-density gas flow which is 
generated as a satellite travels through the upper atmosphere. In partic- 
ular, the use of the proposed Molecular Beam Laboratory will be considered 
for providing the required beam definition and orientation. Engineering 
and scientific gains would be generated by the results of this experiment 
which utilizes an aspect of the orbital flight environment not easily re- 
produced in ground based facilities. The study would evaluate means for 
measuring the aerodynamic forces on a selection of surfaces having a broad 
range of material and physical properties. The study would determine the 
desirable number of surface materials, the range of surface temperatures, 
the range in degree of surface contamination, the number of surface coatings, 
and the angles of attack to be tested in the proposed experiment. Finally, 
the feasibility would be determined of correlating the force measurements 
with changes in gas properties. 

Justification 

The determination of the feasibility of the proposed experiment is 
desirable in view of the potential benefits the experiment would provide. 


100 



At present, satellite aerodynamic properties cannot be predicted accurately 
because of the lack of information the proposed experiment could provide. 
Analytic studies reveal that satellite aerodynamic properties are a strong 
function of the character of the force caused by the gas surface inter- 
action. The gas surface interaction, in turn, is expected to be a strong 
function of surface properties and surface orientation. The design of 
satellites to take advantage of (or to reduce) the aerodynamic forces and 
torques has not been possible because of the lack of information on the 
forces caused by the gas surface interaction. Knowledge of the character 
of the surface forces and the major influences on these forces is necessary 
to the design of satellites to have specified aerodynamic drag, lift, and 
torque characteristics. Such knowledge is also needed in order to inter- 
pret the dynamic response of satellites in the atmosphere as is done in 
the determination of atmospheric density from satellite drag measurements. 

In addition to the engineering information provided by the proposed 
experiment, the results would also contribute to the scientific knowledge 
of the gas surface interaction at satellite velocities. The expected 
results would compliment both orbital and ground-based molecular beam 
studies which provide force information indirectly. The proposed experi- 
ment would then serve to guide future investigations into the more subtle 
details of the interaction. 

Finally, the study of the feasibility of the proposed experiment 
is justified on the basis that the experiment may be a relatively in- 
expensive method of obtaining valuable information. Consequently, the 
experiment could require few equipment components with low development 
cost and short development time. The information gained would be of 
immediate engineering value and would be a valuable complement for future 
gas-surface experiments. 
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Method 


The objective of this study is to determine the feasibility of 
performing a space experiment to measure the aerodynamic forces on sur- 
faces as a function of gas and surface properties and surface orientation. 
The study is divided into three areas: (a) Methods of making measurements, 
(b) Selection of surfaces and surface conditions, and (c) Correlation of 
results with gas properties. 

(a) Measurement Techniques 

The feasibility of making the measurements required will be 
Investigated taking into consideration the expected low level of force 
and the perturbing influences of environmental factors. To be considered 
is the feasibility of the aerodynamic forces acting on flat or shaped 
surface samples exposed to the gas flow generated by the motion of the 
satellite through the atmosphere. The perturbing influence of molecules 
reflected from the satellite may require that the surface samples be ex- 
tended on a boom ahead of the vehicle. Methods will be evaluated for 
measuring the forces, orienting the surface samples, and changing the 
surface properties. 

The accuracy of the proposed measurements is to be evaluated 
considering perturbing environmental influences such as upper atmospheric 
wind and density fluctuations. Methods of calibration and monitoring of 
the environment will be considered as means of Increasing the accuracy of 
the proposed experiment. 

(b) Selection of Surfaces 

The selection of surfaces to be tested in the proposed ex- 
periment will consider the need to reduce satellite payload weight and 
volume while yielding results of the widest possible interest. The se- 
lect! cn of surfaces will be on the basis of providing information of the 
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many factors thought to Influence the forces caused by the gas surface 
interaction. Among the factors of interest are the influence of surface 
material, surface temperature, surface roughness, surface coatings, sur- 
face contaminates, and surface angle-of-attack to the flow. The feasibility 
study would establish a series of experiments which best Isolate the in- 
fluence of the individual factors. The surfaces selected will span those 
used in satellite construction so as to provide engineering Information 
for future design. 

(c) Measurement of the Influence of Gas Properties 

Since the gas surface interaction is Influenced by both the 
surface properties and the gas properties, the feasibility is considered 
of determining the influence of the gas properties on the surface forces. 

The upper atmospheric gas composition, temperature, and degree of ioniza- 
tion are a strong function of altitude, geocentric latitude and longitude, 
and time. To be investigated is the possible correlation of the measured 
forces with the changes in gas properties that occur naturally over the 
orbit. The feasibility will be studied of identifying the gas-property 
Influences on the forces caused by the gas surface interaction. Study 
will be made of the orbital parameters which provide the best conditions 
for the experiment. 

Personnel 

The principal investigator of the proposed study is Dr. G. R. Karr 
(resume attached) who has done considerable work on the theory of the gas 
surfacr interaction and satellite aerodynamics. Since Dr. Karr has pri- 
marily theoretical experience and capability, there is a recognized -peed for 
cooperation with personnel who have experimental capability and experience. 
In view of the good working relationship which exists between Dr. Karr, 

The University of Alabama in Huntsville, and NASA Marshall Space Flight 
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Center, Huntsville, it is proposed that the theoretical expertise of 
Dr. Karr be complemented with the experimental expertise of MSFC personnel 
such as Dr. P. Peters (a surface physicist) and/or Dr. R. Smith (an atmo> 
spheric physicist) both of the Space Science Laboratory at MSFC. 
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APPENDIX 


Listing of computer programs developed and employed under NASA Contract 
Number NAS8-28248. 


Program Name 

Page 

LESQA 

106 

RHORAT 

111 

AFILIP-HIGH CM 

119 

AFILIP-LOW CM 

123 

CDCLEV 

127 

CFEVAL 

130 

RUFSPH 

133 

CLELAN 

137 
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PROGRAM LESQA 


This is a computer program which analyzes denf<lty, temperature, 
velocity and altitude measurements from falling sphere experiments. 
Also in the input data are the densities and drag coefficients em- 
ployed by the original experimenters so that acceleration data can 
also be deduced. Both ascent and descent data are eai>loyed. The data 
are employed in an orthogonal curve fitting routine so that ascent and 
descent data can be correlated at equivalent altitudes. Speed ratio 
effects, molecular weight' changes and drag coefficient variations are 
all taken into consideration. A density Is determined which best 
represents the data based on the measured properties and those pre- 
dicted by theory. A density profile is thus determined using both 
ascent and descent information. 
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«KUN»/ IP LtSQA»l)AHXXX)(XXXXX»tNVlP-ljYN-(jK»b» IbO 


XXMb)= .399+.ulb*S 

GSlFtS) = .997b6l7977-.0l842b96f.*S 

DIMENSION C3D(9nu) *C3A(900) »HHO(9nO) »C0CD<9n0) »nACHO(QOO) »C0 Fmd(90 
10) »SK0(9U0) »CDCA(900) » DACHA (900) »CDPMA(9U0) »SRA (900) »VrATIO( 9O0) 
DIMENSION VAORtqOO) »VDOP(900) 

DIMENSION WVA(90U) »WVD(900) >CVA(90U) »CVD(900) *ALPHVA(900) >BETVA(90 
10) »BETVD(900) »ALPMvD(900) 

DIMENSION AT(900) » TEMP (900) »k>T( 900) *CT (900) • ALPHA! (900 ) # beta T( 900) 
1»TOR(900) 


dimension mD(9U0) 

DIMENSION ALT (SOU) » YAOK ( 90(J ) » YPOR ( 900 ) 

DIMENSION taA(900) »CA(900) *ALPHAA(900) »BETAA(900) » T1 ( 900 > * T2 ( 900 ) t 
173(900) »YA (900) » YD (900) »CD(900) *ALPHAn(900) *BETAD( 900) »RD(900) »PA 
1 (900) 

DIMENSION AA(900) »KH0A(900) »CDA(900) *RHOO(900) » CDD ( 900 ) » AD ( 900 ) * 
1TEMPA(900) »TEMpr)(900) »VA( 900 ) »Vn(900) 

10 READ (5002) (SA) 

302 FORMAT (13) 

REaD(5»250) (C3»E) 

REmD(5»301) (MA*MD»N»ALTM) 

READ (5» 300) ( AA ( I ) » rHOA ( I ) , TFMPA ( I ) »CDA ( I ) , vA ( I ) » Ir l » mA ) 

READ (5» 300) (AD(T) »hHOD ( I ) » TEMPD ( 1 ) »CDD(1) » VD ( I ) » 1=1 »M0) 
WRiTE(6>251) (C3»E) 

V«RlTE(fa»304)MA»MD»N» ALTM 

WRITE (b 003) (AA(1) fRHOA(l) *TEMPA(I) »COA(I) »VA(I) ririfMA) 

WRITE (6 003) (AD( 1) »RH0D( 1) »TEMPD( I ) » CDD ( I ) » VD ( I ) » 1=1 »MD> 

304 F0RMAT(///3HMA=» 14»3HMD=» I40H N= » 14 * SHAlTM= »E 1 0 . 5 ) 

303 FORMAT (4X »Ftt,4»2X»Fe.4»2X»F8.4»2X»F8.4*2X»F9.4/) 

6AMMA=1,4 

PI=3. 141592653 
C0NVSM=((jAMMA*.6)**,5 
DO 80 I=1»MA 
AT(I)=AA(I) 

60 TEMP(I)=TEMPa(I) 

DO 81 1=1»MD 
AT (MA4I )=AD( 1 ) 
ei TEMP(MA+ 1 )=TEMPD( 1 ) 


VT=MA+MD 

300 FORMAT (F7.3»Efa.3»F6.4»F5.3»F8.?) 

301 F0RMAT(313»Fb.2) 

DO 70 I=1»MA 

70 WA(I)=1.0 

DO 71 I=1»MD 

71 WD(I)=1.0 
WRITE (6*305) SA 

305 FORMAT (4X* 17hS0DNDlNb NUMPER =*I3/) 
WRlTE(b*201) 

DO 2 I=1*MD 
YID=HHOD( I ) *CnL( 1 ) 

SXD=SXD+AD( I ) 

WR1TE(6*100)AD(I) *YID 
Y = ALOG (HhOD ( 1 ) *CCD ( I ) ) 

YD(1)=Y 
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2 CONTlNUt 

twRlTt <b»c00) 

DO 1 I=1»MA 
YlA=KHOAiI)«CDA(i) 

Y=ALCG(RH0A(I)*CDA(I) ) 

YA(I)=Y 

teRlTL(6»lOO)AA (I) » YI A 
1 CONTINUE 

WRITE (6»2U0) 

DO GO 1=1»MA 

60 WRITE (6»IU0) AAU)»YA(I) 

WRITE (6»2U1) 

DO 61 l=l»MO 

61 WRITE (6»100) ADiUrYOd) 

N=bl 

DO 30 K=1»N 

30 ALnK)=ALTM4.5*(K-l)-»lO.O 
LT = 0 
JT=0 
KF1=3 
KT=3 
KD=2 
KA=2 
KFA=2 
KFU=2 
KCA=2 
KCU=2 
LA = 0 
LD=0 
JA = 0 
JD=0 
KVA=3 
LVA=0 
LVU=0 
KVU=3 
JVA=0 
JVU=U 
KFVA=3 
KFVD=3 

CALL ORTHLb ( AT » TtMP » WT »MT *LT » JT »CT » ALPHAT »BETAT »KT » T 1 • T2 , T3 » INOlT 
1 ) 

CALL ORTHLb ( A A , Y A • W A » MA » LA » JA #CA » ALPHA A » BFT A A » K A » T 1 » T2 » T3 » INUI A ) 
CALL ORTHLb ( AD » YU » wD » PU »LO » JD » CD» ALPHAU»BFTAD »KD » Tl » T2 » T3 » INDID ) 
CALL ORTHLb ( AA . V A r WVA »MA »LVA , JVA »CVA , ALPHVA » BETVA »K VA » T1 » T2 » T3 » IN 
IDIVA) 

CALL ORTHLS(AD»VU»WVU»MD»LVn» JVn»CVD»ALPHVD»BETVD»KVD»Tl»T2»T3» IND 
IIVU) 

11=0 

WRITE (6 » 103) IND1VA» II»CVA( 1) » (I I » C VA ( 1 14 1 ) , ALPHV A ( 1 1 ) * BETVA ( 1 1 ) » 
111=1 »KVA) 

11 = 0 

WRITE (6»lu3)lNDly/U»lI»CvD(l) » (II»CVD<II41) »ALPHVD(ID »BETVD(I1) » 
1II=1*KVD) 

11=0 

WRITE(b»103) 1NU1A»1T»CA<1) » (1I»CA(1I41) • ALPHA A ( I 1 ) » BETaA ( I I ) »TI=1» 
IKA) 

11 = 0 


108 



MHlTk(6»103)IN01U»iI»CU(l)»(XI*C0(ll4l), aLPHAO ( T I ) t BET AO ( 1 1 ) * 1 1 = 1 » 
IKD) 

11=0 

WRlTt<6»103) INUlT»iI»CT(i) » (1I»CT(1I^1) » ALPHAT ( T I ) »BETAT ( 1 1) »Tl = i 
1»KT) 

103 F0KMAT(//13///2X» 1h» 2X»F20.7»/(?X» l<4»2X»3E2n.7> ) 

call COEFS (jA#CA»mLPHAA,BETAA»KCA»BA»T1»T2»T3» IND2A) 
call COEFS (jn»CO»ALPHAn»BtTAn»KCD»PD»Tl »T2»13* IND2D) 

CALL FITr (ALT »N»wA»CA»ALPHAA»BETAA,KFA» YAOR fTl »T2» 1N03A) 

CALL FITY(ALT»N»oL»CD»ALPHAD»RETAD»KFD»YD0R*T1»T2# IND3D) 

CALL FITY (ALT »N»uT»CTf ALPHAT »BETAT»KFT» TOR »T1*T?» IND3T) 

CALL FITY(ALT»N» JVAfCVAf ALPHVA»PETVA»KFVA»VA0R»Tl»T2r IND3VA) 

CALL FITY ( ALT •N»JVU»CVO»ALPHVn»PtTVD»KFVO»vOOP»Tl»T2r lNn3VD> 

^RlTt (6r5U0) (K » ALT (K ) » TOR ( K ) » VAOK ( K ) » VUOR (K ) » K=1 » N ) 

500 FORMAT (//3X»1HK»5X»6HALT(K) *dX»6HT0R(K) » ttX »7HV AOR ( K ) dX » 7HVD0K ( K ) »/ 
1 (2X» 13»Fl0.4»3Ll5.tj) ) 

WRITE (6*450) (I»bA(I) *I=1*KCA) 

WRITE (6*4bl) lT*bU(I) *I=l*KCn) 

450 FORMAT ( //4X » IFi 1 » IbX * 5HPA ( I ) » / ( 2X » 1 5 * 3X * E20 . 7 ) ) 

451 FORMAT (//4X*lHl»lbX»5HBu(l) »/(PX»I5*3X*E20.7) ) 

WR1TE(6*202) 

DO 3 I=1»N 

EANM=24.b8'»’.1235*ALT( I )-.000875»ALT ( 1 )*ALT(I ) 

RG=fl3l4.34/EANM 

SRD(1)=VU0R( I)/(2.*R6*10R(i) )**.5 
IF (SPD(l)-1.2b)lSb»19b*l9b 

195 S=bRU(I) 

CMU = S/LONVSM 
GSl = GSIF(CMD) 

CDF MO ( I ) =(2./(P1**.5) )*(8./(3.*S)+P.*S/15.-8.*S*5*S/2l0. )#GbT 
GO TO 197 

196 S=bRU(i) 

CMU = S/LONVSM 
GSl = GSIF(CMD) 

SI2=l./(b*S) 

Sl4=bI2*bIii 

CDFMD( 1 )=2.^(1 .+S12-.Pb*Sl4)*GSI 

197 DALHD( I )=SRD(I )/CunVSM 
CM=DACHD( I ) 

CDLA( I )=.92^. lb6U7l4/DM-.3b607l4/(OM*nMADM) 
VRAT10(I)=VAOR(I)/VPOR(I) 

SRA( 1 )=VAOR( I)/(2.*RG*T0k( 1 ) )**.5 
IF (SRA(l)-l.Pb) 95*95*96 
95 S=5RA(I) 

CMA = S/LONVbM 
GSl = GSIFICMA) 

CDFMA( I ) =(2./(Pi»*.5)) • (8./(3.*b)48.*S/l5.-8.*S> i*S/2lO . ) *GSI 
GO TO 97 
9b S=SRA(I) 

CMA = S/LONVSM 
GSl = GSIF(CMA) 

SI2=1./(S*S) 

Sl4=SI2*Sld 

CDFMAd )=2.*(l.+S12-.2b*SI4)*GSI 
97 DALHA( I )=SRA( I )/LO'mV.'.vi 
nM=DACHA ( I ) 

COLD ( 1 ) = .92+. IbfU /l4/DM-.3b607l4/ (UM*nM*DM) 
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PHOCUA=EXP(YAOK( 1) ) 

RHOCDD=EXPlYOOR ( 1 > > 

DENO^*=CDt A 1 1 ) < I )-CDCU ( I ) *CUFMA ( I ) 

PHO(l) = (RHOCDA*(CUhyU(I>-CDCDn > )-RHOCDD*(cnFyA(I)-CDCA<I) ) )/nENOM 
RH0IN=1./RH0(I) 

XX= (RH0CDA*RH01N-tDCA ( I ) ) / (CCFPA ( I ) -COCA ( I ) ) 

IF (XX)77»79»79 

77 C3A(I)=-1. 

GO TO 78 

79 X = XXF(CPA) 

C3A( I)=-(RM01N**X )*AL06(XX) 

78 XY=(KH0CUD*RH01N-CUCC(1) )/(CDF^Pn)-CnCD(I) ) 

IF(XY) 87»69»89 

87 C3D(1)=-1. 

60 TO 88 

89 X = XXF(CPb) 

C3U( 1 )=-<RHOIN**X )*ALOG(XY) 

88 AOVPD=RHOCUA/RhOtUU 
FPKA=CUFPA 1 1 )/CnFMb< 1 ) 
tONRArCDcA 1 1 ) /CDCb 1 1 ) 

TR'IN=( ( l.-AOVRb)/tAOVRb/lVUOR( I )-l./(VAOR( !)♦*?. ) ) >/(iJ.*Pb) 

V*R1TE(6» 102) ALT ( 1 ) » Y AOR ( 1 ) » RFOCHA , YPOR ( I ) » RHOCDO » FMRA » AOVRD » CONP A , 
ITPIN 

3 CONTINUE 

WRITE (6» 333) (ALT(i) »SRA(1) »CCCA ( I ) »CUFMA ( I ) »DACHA( I ) »I=l»N) 

WRITE (6* 334) (ALT (1 ) *SRb( I ) »CDCD ( I ) *COFPD (I ) »DACMD( I> » 1 = 1 »N) 

WRITE (6 » 33b) (ALT U ) »C3A( 1 ) »C30( I ) » VRATIO( I ) »RHO( 1 ) , 1 = 1 »N) 

333 F0RPAT(//4X»3HALT»9X»3HSRA»qx»4HCDCA»9X»5HCDFMA,8X»5H0ACHAr/(pX»Fl 
10.4.4E13.5) ) 

334 FORMAT(//4X»3HALT»10X»3HSRD» 10X»4HCDCD»9X»5HCnFMb»8X»5HDACHD»/(?X» 
1F1 u.4»4E13.5) ) 

335 format (// 4X»3HALT»10X»3HC3A»inX»3HC3U»8X»6HVRATT0»8X»3HRH0»/(?X»F I 
10.b*4E13.5) ) 

60 TO 10 

100 FORMAT(2X*F10.4f 2X»E12.6) 

102 FORMAT (2X»F10.4»titl4.b) 

200 FORMAT (///bX»5HAA( 1 ) » 10X» 3HYIA//) 

201 F0RMAT(///bX#5HAb(l ) » lUXr3HYin// ) 

202 FORMAT (///bX»3 H alt » 12X » 2HY A » bX » bHRHOCDA » lOX »2HYP» lOX » 6HRH0CDU » 8X f 4 
IHFMRA# 10X»bHAOVRD»9X»5HCONRA) 

250 FOHMAT(2E20.10) 

251 FORMAT(//f 4X»3HC3=»F20.10»4X»2HF=»E20.10) 

ENU 
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PROGRAM RHORAT 


This program takes falling sphere data and performs the usual 
analysis to determine the temperature and density profiles. The data 
in the region of overlap of the ascent and descent trajectories is 
treated as if the atmosphere were influenced by isentroplc waves. 

Thus, differences in density at the same altitude will result in 
differences in temperature at that altitude according to the isentroplc 
relations. 

The program also Inputs various standard atmospheres so that 
ratios of the measured density and temperature can be readily com- 
pared to the standard values. 
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«HUN,/TP RhOKAT,UAHXXXXXXXXX,ENVIH-0YN-6K»5,l5U 


XXF(S)“ .199^.016*S 

■ . 99756 1 7977-. 018*«2<.966*S 
tANHF(G>“2‘<.6B^.l235«6-.00087M«(;«G 
parameter NPP-200 

OlMERblON M(5U0I * NAD (500) iNAA(^OQ) ,TU(bOO) »TU2(S00 ) .TmeISOO) i 
1 SiG( 500) ,PU(500I ,PU2(S00) ,PME)500) .qjPtbOO) iTmSRCsOO) ,SiSR(500) 
DIMENSION RHDD( 1 00 ) . R A I SO ( 1 00 ) , T R A ( ] 00 ) , C I SO ( 1 00 ) 
dimension TTA(NPP|,»«TA(NPP),CTA(NPP),APTTA(NPP),BETTA(nPP) 
DIMENSION GR(lOO)t<>^A(100)tALPHG(IOn)tBETA(i(IOO),6RD(IOO)iC6(100) 
DIMENSION TINA(NPP),TINO(NPP)iThA(NpP),THO(NPP) 

DIMENSION SRA(NPP),COFMA(NPP),OACHA(NPP),CDCA(NPP),SRO(NPP),COFMD 
1 (NPP),0ACH0(NPP) ,CDC0(NPP) 

DIMENSION AA(NPP),RHOA(NPP),COA(NPP) ,RMOD(NPP) ,CDD(NPPJ,A0(NPPJ, 
1TEMPA(NPP|,TEMPD(NPP),VA(NPP),VD(NPP), 

IRH0GA(|00),NAN(I00)|RH0RA(100), COGA(100)|CORA(IOO)*RhOGD(IOO)|NON 
1 ( IOOI.RHOROIIOO) ,COGD(lOO)iCOROriOO) 
dimension RHL5(NPP),RHLSA(NPP),RHLSd(NPP) 
dimension RHOSA (NPP) .RHOSO(NPP) 

DIMENSION ALS(NPP),RH0S(NPP),TES(NPP),TI(NPP),T2(NPP),T3(NPP),FP8( 
INPP),SM0L(NPP|iW(NPP),CS(NPP),AlPS(NPP)i8ETS(NPP) 

DIMENSION AT(NPP) 

1 ,XA(NPP) ,XO(NPP) ,C3A(NPP) ,C30(NPP) 

1,HHA(NPP),CRA(NPP),APA(NPP),BEA(NPP),RHD(NPP),TEMAD(I00) 

1 .TTAINPP) ,TTD(NPP) 

CC ■ 307O72C00. • . 6H 

kta ■ a 

JTA ■ 0 
LTA ■ 0 

HEAD (5,601 ) (MS) 

READ (5,800) ( ALS( I ) ,RHOS( I ) ,TES( I ) ,FPS( I ) ,SM0L( I ) , I«I ,MS) 
head (5,600) ( GR ( I ) , 1 ■ I ,MS ) 

wR I TE ( 6,599 ) 

599 F0RMAT(IH0 ,8X , • alt • ,8X , *GRAV I TY • ) 
mRITE (6,601) ( ALS( I ) ,GR( I ) , I>I ,MS) 

600 format (I0X,FI0.5) 

WRITE (6,803) (MS) 

WRITE ( 6 ,802 ) ( ALS( I ) ,RhOS( I ) , TES ( I ) , FPS ( I ) ,SM0L ( 1 ) , I-I ,MS ) 

flo 1 format ( I I 0 ) 

800 FORMAT (FI5.5,EI5.5,FI5.5,EI5.5,FI5.5) 

803 format ( /// ,2X , 3HMS- , I 10,2X , I 3HU.S. ST AND A RD , / / , 9 X , 3h A lS , I 2 X , MHRhO 
1S,12X,3HTES,I2X, 3HFPS, I 2X , HHSMOL i // ) 

802 FORMAT (2X,FI5.5,'>*I5.5,FI5.5,E15.S,f 15.5) 
tR'.OOOI 
DO 60 I ■ 1 , MS 

60 HHLS( I ) «AL0G(RH0S( I ) ) 

KKK>6 

call ORTmLS (ALS,RHLS,W,MS,O,0,CS,AlPS,BETS,KKK,TI ,T2,T3, INDl ) 

GAMMA> 1 • H 

C0NVSM«(GAMMA«.5)**.5 
PI*3. IMI592653 
DO 8M I -80 , I 20 
M( I )«0 
TUI 1 ) >0. 

TU2( I ) -0.0 
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TML U ) *0.0 
si<j( I ) >u. 

P J ( 1 ) >0 . 

< I) *n. 

PML ( I ) *n . 

8 M s I P ( I » ■ 0 . 

10 KL AU ( b . 302 ) ( SA ) 

HL AO ( b . 3 1 0 ) ( MAU I 
WtAO (S,20n> (C3.E) 

WLAO(S,301 )(HA|M0,N,ALTM) 

Wf. A0(S,3 00> I AA( I ) ,HhOA ( I ) ,TfcMPA( I ) ,cOA( I ) ,VA( 1 ) , 1*1 ,Ma » 
Nl.AO(b,300)(A0(I),HH0D(n,TEHP0(I),cn0(l),VD(l),l»l,M0) 

tea ■ tempa n ) 

TED ■ TEMPO ( 1 ) 

00 MR I » I ,MA 
NAA(I|«AA(I)>*2 
ATII) « AA(I) 

00 8S I * 1 ,MD 
TEMAO ( I ) -OtO 
R H 0 0 ( I ) ■ 0 • 0 
8b NAJ(I)3A0(I)^.2 
00 M7 1 « 1 ,b 

•49 wTA C 1 ) ■ .b 

00 M? I » 6,MA 
•47 A ( 1 ) = I .0 

CALL OHTHLS 1 ALS ,6R , <4 . ms ,0 ,0 ,CG , ALPhG i8ET aG . *4 , T I , T2 , T3 , 1 no I ) 
CALL pity ( a a ,MA ,0 ,CG , ALPHG ,BET AG , *4 ,6«A , T 1 , T2 , I N03 1 
CALL FITY ( AO ,M0 ,0 ,CG , ALPHG ,RE T AG , *4 ,GRD , T I , T2 , I N03 ) 

601 FOmMAT (RX.ZFlU.b) 

WRITE (6.6011 ( AA ( 1 ) ,GRA ( i ) , I > 1 ,MA ) 

mRITE (6,601) ( AD ( I ) ,GRD( I ) I I > 1 ,M0 ) 

00 b I « I • 1 00 
COMA ( I ) «0.0 
b CORO ( 1 ) *0.0 

WRITE ( 6 . 30b ) bA 
WR 1 TE ( 6 • 3 1 1 ) ( MAD ) 

WRITE (6,201) ( C3 , t ) 

«RITF(6,30‘4)MA,M0,N,ALTM 

#»RITF(6,303)(AA(I),Rh0A(1),TEMPA(I).CDA(I),VA(I),I«1,hA) 
wRlTE(6,3n3)(M0( 1 ),RHOO( I ) ,TEMPD( I ) ,C00( I ),VD( I ),I*1 ,MD) 

30b format (*4X, 17HS0UN0ING NUMBER »,I3/) 

301 F0RMAT(3I3,F6.2) 

300 FORMAT ( F 7 . 3 I E8 . *4 ,Fb . *4 ,Fb. 3 , F8 . 2 ) 

311 FORMAT ( *4X , *4HMAD* , I *4/ ) 

3 1 0 FORMAT (12) 

30*4 F ORM A T ( ///3hMA ■ , I *4 , 3hmO* , 1 *4 , 3h N ■ , M , 5h A L T M ■ , E 1 0 . b ) 

303 FORMAT ( R X , F 6 . 2 , 2 X , E 7 . 2 , 2 X , F *4 . Q , 2 X , F 5 . 3 , 2 X , F 6 . 1 / ) 

302 format (13) 

200 FOMMA T ( 2E20. I 0 ) 

201 format (//, *4 X , 3HC 3* ,E 20. 1 0 , *4 X , 2hE », e 20. 1 0 ) 

call FITY ( A A ,MA ,n ,CS , alps ,RF. TS ,X<K .RHLSA , T 1 , T2 , 1 ND3 ) 
call FITY ( ad ,M0 ,0 ,CS , alps ,RF TS ,xKX .RHLSD . T 1 , T2 , 1 N0*4 ) 

MNU*Q 

00 bb 1 ■ 1 ,M0 
bb RHOGO ( I ) »RhOO ( I ) 

76 00 1 2 I ■ 1 , mA 
AL T A « A A ( 1 1 
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R(i«831H.3M/LANMF(aLTA) 

SRA< I )»VA< I )/(2t»RG*TEMPA( I I )*«,5 
IF (SRA( I 1-1 ,2bl 9b.9bi96 
9b S»bRA ( I ) 

CMa ■ 5/CONVSM 
Gb i ■ gS I F ( Cm A I 

CUFmA( I ) ■(2,/(PI**,S) )«(gt/(3,*S)4g,»S/lb,.0,«S»S»S/21O, )«GSI 
GO TO 97 

96 S»5R A ( I » 

CMA ■ S/CONVSM 
GbI - G5IF (CMA I 
S12«I ./(S*SI 
blM«5l2«SI2 

CUFMA( I )-2.*( l.*Sl2-.2b*SI9)*GSI 

97 dacha ( I ) »SR A ( I) /CONVSM 
DM-OACH A I I ) 

CDCA( 1 )• » 92 * • 16A071S/DM-#36A07IM/|Um«OM*DM) 

OC-COFMA ( I » -CDC A ( I ) 

NN»20 

RHOO-HHOA ( I ) 

K»- 1 

XA ( 1) ■ XXF (CMA ) 

X » XXF ( CMA ) 

C3A(II • . 2 1 2* ( CC/t ANMF ( alt A ) ) • »x A ( n 

C3 ■ C3A ( I ) 

30 F X-RHOA ( I I »CDA ( 1 ) / ( COCA ( II ♦DC^EXP ( -C3« ( RHOO»*X HI 
CALL wEG I T ( RHOO ,FX ,t ,K ,NN I 

GO TO ( 30,3 1 , 32,33) ,K 

3 1 RHUGA ( I ) >RHOO 
GO TO 3M 

32 RHOGA ( I I «0.0 
GO TO 3M 

33 RHOGA ( I ) ■- 1 . 0 
3M HAN ( 1 ) »NN 

RHORA( I )»RHOGA( I )/RHOA( I I 

COGA( I )«COCA( I )4.DC»EXP(-C3*(RH0GA( I )**X ) ) 

KHOSA ( I I »EXP ( RHLSA ( I I I 
CORA ( 1 I-RhOGA ( I )/KHOSA( I I 

12 continue 

I F ( M ad-50 I 20,21,21 

20 00 58 JX > 1 , M AO 

TfcMAD(JK)«TEMPA(MA-MAD^JX) 

RHDD(JK)«RH0GA(MA-MA04JK) 

58 TEMPD(JK)»TEMAO(JK)«(RHOGD(JX)/RHDO(JK))»«(tR) 

GO TO 22 

2 1 TLMPD ( I I ■ TEO 

22 00 II I ■ 1 ,MD 
ALTO-AO ( I I 

RG»831R.3R/EANMF(ALT0) 
bHO( I )«V0( I )/(2 .»RG»TEmP 0( I I (••♦5 
IF (SRD( I l-l .25) 19S, I9b, 196 
I 9b S-bRO ( I ) 

CMC » S/CONVSM 
GbI « 6SIF ( CMO ) 

C0FMD( I ) »(2./IPI««.b) )»(&./(3.»S)^ft.«S/lb.-8.«S*S*S/2lO. )«GSI 
GO TO 197 
I 96 S»bRO ( I ) 
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c«0 • S/cONVSM 
Gbl ■ gSIF (CMU ) 

SI2-1 ./«S*SI 
5IH»S12*SI2 

COfMO( 1 )-2.*« l.*S12-#2b»SIH»*G5I 
197 DAChO ( I » -SRO U I /CONVSM 
OrisUACHO ( I 1 

COCD(I)«.92^.1660;M/OM-,36A0714/(Om*DM«OM) 

OC«COFMD ( I 1 -COCO ( I ) 

NNb2Q 

K«- 1 

RHOGbHHOO ( I 1 
XD ( I ) - XXF ( CMD ) 

X ■ XXF 1 CMD ) 

C3D(I) » .2l2*(CC/e.ANMF( ALTO) )**XD( I ) 

C3 « C30 n ) 

MO F X «RMOO ( I ) ‘COD ( 1 ) / ( COCO ( 1 ) ♦DC^E XP 1 -c 3» < f’HOO** X ))) 

call weg I n rhou X ,e ,k ,nn ) 

GO TO 1 MO I H I , M2 I M3 ) ,K 
M 1 RHOGD ( I ) *RHOO 
GO TO MM 
M2 HMOGD I I ) *0.0 
GO TO MM 

M3 RHOGOl I) ■- I .0 
MM NON 1 1 ) »NN 

HHURO( I )>RH0G0( I )/RHOD( I ) 

C0GD(I)-CDCD(1)^DC*EXP(-C3«(RH0GD(I)»»X )) 

RHOSO( 1 ) »EXP(KHLSD( I ) ) 

COHO( I )>RHOGO( I )/RHOSD( I ) 

RAISOI I )»JRH0D1 I )/HH0G01 I ) )**.M 

1 I CONT I NUE 

WHITE (6.902) 

white (6,900) ( A A ( I ) .HhOGA ( I ) .RhOHA ( I ) • C 0 G A ( I ) , C 0 H A ( f ) , N A N ( I ) , 

1 X A ( I ) ,C3A ( I ) , 1 ■ I .MA ) 

WHITE (6.901 ) 

WHITE (6,900) ( AD( I ) ,RhOgO( I ) ,RhOO( 1 ) , CDgO ( I ) , CDhO ( I ) , NON ( I ) , 

1 X0( I ) ,C30( I ) . 1*1 ,M0) 

RB*0*0 

DO 61 X-2.MA 

8B»BR-(AA(K)-AA(X-i ) )»(RHOGA(K)*GRA(K)-HHOGA(K-i )*GHA(X-1 ) l/ALOG 
1 (HH0 6A(K )*GHA(K)/(HH0GA(K-1 )*6 Ra(K-1 ) ) ) 

61 T I NA ( K ) -OB 
BBsO.O 

DO 63 K>2,MD 

BB»BB-(AD()(:)-A0(K-1 ) )»(RH0G0(K)*GHD()C1-RH0GD(K-1 )»6RD(K-1 ) 1/ALOG 
?(HhOGO(K)«GRO(K)/(RHOGD(K-1 )*GHD(K-I ) ) ) 

63 T I NO ( K ) *BB 

THA ( 1 ) «TEMPA ( 1 ) 

DO 6M KK«2,MA 
alt A aAA ( KK ) 

RG*83IM,3M/EANMF(ALTA) 

6M THA(kk)*TINA(KK)*10U0./(HH0GA(kk)»R6)*RH0GA( 1 )«THA( | )/rHOGA(kk) 

THD ( 1 ) -TEMPO ( 1 ) 

DO 6B XX-2.MD 
ALTD-AO ( KK ) 

H««831M,3M/EANMF(ALT0) 

6S TH0(FK)«TIND(KX)*1000. /(HMOGD (KX)«Rs)*HHOGD( * )*THD( 1 )/RH0G0(fK) 
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DO b7 I • 1 , mO 

THA( I )»TLMA0( I »/TEMPD( I ) 

57 C1S0( I |.HA1S0( I )/TRA( I ) 

write I6,610» «AA( I > .TEMPa« I ) tlHA( n .!■» iMa» 
WKITE(6,611MA0(i).TEmPD(i)|Th0(i)|RAiS0(i)»TEmAD(i),CiS0(i 
J ) , I ■ 1 ,HD ) 

00 71 I ■ 1 , HA 

01pA-ABSI(TEmPA(i»-Th*<I»»/ThA(i))^0IF^ 

7 1 TEMP A I I ) - ABS ( Thai I ) ) 

IF (HaO-50) 23,2M,2‘4 
2J HA1»Ma0^1 
GO TO 25 
2m MAI > 1 

25 00 72 I -MA 1 ,M0 

DIFD»AHS( ( TEMPO! I )-TMD( I ) )/THQ( 1 ) J^OlFO 
72 TEMPO! I J»ABS« ThD( I ) ) 

SUMO-0 I F A * FO 
TEMPa 1 1 ) ■ TEA 
IF (SUMO-ER) 7M,7«!,75 
75 MNO»MNO^i 

IF lMNO-5) 78,78,7M 
?a GO TO 76 

7h continue 

00 HI 1 - 1 t MD 
JC«NAD ( I ) 

M IJC » «M ! JC ) ♦ I 

TU2IJCI»TU2(JC)^TEMP0( I )«*2 
PU2(JC)»PU2(JC)^C0R0! I )»«2 
TU( JC)-TUIJC)^TEMPD! I ) 

PU ! JC I «PU ( JC ) ♦CORO ( I ) 

81 continue 

DO 62 I - 1 ,MA 
JC-NA A ( I I 
M ( JC I »M ( jC ) ♦ 1 

PU2(JC)»PU2(JC)*C0RA( 1 )«*2 
TU ( JC) »TU ( JC ) ♦TEMPA ( I ) 

PU ! JC ) -PU { JC » ♦CORA I I I 
TU2(JC)»TU2{JC)^TEMPA(I)«»2 

82 CONTINUE 

00 83 I >80 . 1 20 
TML ! 1 ) -TU I I ) /M ! 1 ) 

Sl&l I )■! ! TU2( I I » I-TME! I )*»2)««.5 
PMt ! I ) «PU ( I I /M ( I ) 

J ■ 1-79 

TMSR ! 1 I « TME ( I ) /TES ( J ) 

SI5R!I) ■ SIG(1)/TES(J) 

«3 SIP1I)-1(PU2(1)/M(I))-PME(I)**2)*«,5 

WR1TE(A|A5Q)(I.M(I),Tme<1»»5IG<I>iPmE(I)i5IP(I».TmSR(I),SISRIi)» 

II -80.1201 

650 FORMAT ( 6X . 1 H I i6X I MHM ( I ) . 7X . 6HTME ( I ) . 9X.6 hSIG(|)i 9X,6hPHE(I). 

I 9X,6HSIP(I)|9X,7HTM5R(n.9X,7HSISR!I)//(5X.I3|5X,l3,6Eib*7J) 

610 format (6X.5hAAII),10X.8hTEMPA!i).1mx.6HThA(I|//(MX.F10*5.2E2Q*10) 
3 ) 

611 F0KMaT(6X.5HA0<I1»10X,8hTEMP0II)i13Xi8HTema0I1)|9X,8hRaIS0(I>* 
112X,6HTRA!I),1UX,7hC1S0(I)//(MX|F|0.5.5E2U«10)) 

900 FORMAT I 2X .F 1 O . H I me 1 H • 6 . I 1 0 , 2E 1 M • 6 1 

901 FORMAT ( /// I 6X , 3HALT ,9X .5HRM0G0 , 9X »5HRHDD ilOX. MhCDGD . 1 OX , MMCDRO » 
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inx,3HA0N,9X,»XA»,13X,*C3AM 

902 FOHMAT ( /// |6X , 3HALT , 9X .bHRHOGA , 9X , 5HRH0R A , J OX , MHCOQA , I OX , 9HC0R A , 
111X,3HAAN,9X,»X0»,13X,»C3D») 

GO TO 10 
CNO 

UXQT 


4 1 


80* 

1 

.999 

-E5180.65 

9.065 

81 . 

1 

.662 

-E5 1 80.65 

9.888 

82. 

1 

. 382 

-E5 180.65 

5.877 

83. 

1 

. 150 

-E5180.65 

7.067 

89 . 

9 

.563 

•C6 1 80. 65 

8.996 

85. 

7 

.955 

>£6180.65 

1.021 

8 6. 

6 

.6 17 

-£6 1 80.65 

1 . 228 

87 . 

5 

.509 

-E6 1 80. 65 

1.976 

8 8. 

9 

.579 

>£6180.65 

1.779 

89. 

3 

.819 

>£6180.66 

2.133 

90. 

3 

• 1 70 

-£6180.65 

2.563 

91 . 

2 

.598 

-£6183.63 

3.127 

92. 

2 

. 137 

>£6186.62 

3.802 

93. 

1 

.763 

>£6189.59 

9.607 

99. 

1 

. 959 

>£6192.56 

5.566 

95. 

1 

.211 

>£6195.51 

6.702 

96. 

1 

.00b 

198.95 

8.052 

97. 

8 

.915 

>£7201 . 37 

9.693 

98 . 

7 

.099 

>£7209.28 

1.151 

99. 

5 

.9 11 

>£7207. I 6 

1.371 

100. 

9 

.979 

>£7210.02 

1 .629 

lot . 

9 

. 159 

>£7219.86 

1 .996 

102. 

3 

.993 

-£7219.66 

2.316 

103. 

2 

. 995 

-£7229.93 

2.799 

109 . 

2 

.992 

-E7229. 1 8 

3.290 

1 05. 

2 

.117 

>£7233.90 

3.810 

1 OA . 

1 

. 809 

-£7238.58 

9.965 

107. 

1 

.593 

-£7293.23 

5.215 

108. 

1 

. 323 

-E7297.85 

6.07 1 

1 09. 

1 

. 1 39 

>£7252.99 

7.095 

110. 

9 

.829 

>£8257.00 

8.150 

111. 

8 

.360 

>£8266.99 

9.568 

112. 

7 

. 153 

>£8275.85 

1.117 

113. 

6 

. 153 

-£8285.20 

1 .296 

119. 

6 

.321 

-£8299 . 92 

1 .996 

1 IS. 

9 

.623 

>£8303.78 

1.719 

116. 

9 

.035 

>£8313.01 

1 .966 

117. 

3 

.536 

-E8322. 1 9 

2.239 

118. 

3 

>112 

>£8331 .33 

2.590 

119. 

2 

.798 

>£8390.93 

2.870 

1 20. 

2< 

>936 

-E0399. 99 

3.233 

80 . 

9.569 




81 . 

9.56 1 




82. 

9.558 




83. 

9.555 




89 . 

9.652 




86. 

9.550 




8 6. 

9.697 




87. 

9.599 




88 . 

9.59 1 





-E328.96M 
-E328 . 96H 
-E 328,969 
-E32fi»969 
-E328.969 
-E228,96*4 
-E228.96M 
-E 22ft • 96 H 
-E228 . 969 
-F 22« • 969 
-E228.96 
-E228, 96 
-E228.96 
-E228*96 
-F828,95 
-F228.99 
-E828 . 99 
-F828.92 
-F I 28.91 
-F128,90 
-E 128.88 
-E 128.86 
-F 1 28 . 83 
-E128.81 
-E 128.78 
-E128.75 
-F 128.72 
-F 1 28 . 68 
-E 1 28 . 669 
-E 128.60 
-F128.56 
-E 128.51 
28.97 
28.92 
28.37 
28. 32 
28.27 
28.22 
28.17 
28.12 
28.07 
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0UU3H 1 

000 

OUU3H2 

000 

0U03M3 

000 

0U03MM 

000 

nuu3*<b 

000 

0UU3H6 

000 

0U03M7 

000 

0003*48 

000 

0U03M9 

000 

OUU JbO 

000 

0UU3b 1 

000 

0UU3b2 

000 

OUU3S3 

000 

0UU3S9 

000 

0U03^b 

000 

0U03S6 

000 

0UU3S7 

000 

OUU358 

000 

000359 

000 

000360 

000 

00036 1 

000 

000362 

000 

000363 

000 

00036M 

000 

00036b 

000 

000366 

000 

000367 

000 

000368 

000 

000369 

000 

000370 

000 

000371 

000 

000372 

LLU44^ 
V U ii 


89. 

9.538 

9q. 

9.535 

91. 

9.532 

92. 

9.529 

93. 

9.526 

9h . 

9.523 

95. 

9.520 

96. 

9.617 

97. 

9.5 19 

9B. 

9.511 

99. 

9.508 

100. 

9.505 

101 • 

9.502 

102. 

9 . *499 

103. 

9. *496 

I OM . 

9. *493 

105. 

9 . *490 

106. 

9. *488 

107. 

9.M85 

108. 

9. *482 

109. 

9,H79 

1 10. 

9.M76 

111. 

9. *473 

112. 

9. *470 

113. 

9. *467 

119. 

9.969 

115. 

9,961 

116. 

9. *458 

1 17. 

9.955 

118. 

9.952 

1 19. 

9,999 

120-. 

9,997 
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PROGRAM AFILIP-HICH CM 


This program computes values of drag coefficient for various 
values of Knudsen number and Reynolds numbers for Mach numbers above 1. 
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«JHUN,/TPC AFILlP»OAHXXX)(XXXXX,OR8|T,b,HOO 


high C*1 

bFOM.lS MAlNiHAlN 

GSIFJbl • . VV7b6l 7977-.01BH2A96A«S 
XF(S) ■ .399^.0|6*S 

CUfMLF<S>*<2./<RI*«.SII»(8./(3.*S)^A.»S/lb.-8.*S»S*S/210.)»GSI 
COFHMFIS) ■ 2.»l I .♦» ./iS*SI-.25/IS*5»S*Sn»GSI 
CUCF ( 5 » ■ • 92^ t 1 <»607 1 •« /S-. 166071 9 /(S*S«SI 

UIHLNSION XH( lUO) »H£( 100) ,C0( 100) ,B<N( 100) I 
I C0L5( 100) ,CR( lUO) 

PI ■ 3.191b926b 
GAMMA • 1.9 

C0NV5M « I gamma* . 5 ) •• *b 

10 HLAOf bi 1 *eND«lUO)CH,N| 

1 FORMATIF lOtOi lb) 

R£A0(b|33) XLiCLtGSl 

33 FORMAT! 3F IStO) 

RLA0(b.2) (XM(I),RE(I)»CD(I),t>l|Nl) 

2 format ( 3F 20*0 ) 

GbI > GSIF(CM) 

SLLK ■ n.o 
SLL<2 ■ 0.0 

ST • n.o 

SYKN . 0.0 
WH I TE ! 6 , 1 1 ) 

11 FORMAT! IHI ,13X, ‘PRESENT EXPERIMENTAL SPMERE OATA‘1 
NR I TE ! 6 ,S )CH 

5 format ! IHO,/ ,/ .25X , »CM * * . F 6 , 9 , / , / , 1 I X , • X M • , I 7 X , ♦ RE • , 1 9 X , • CO » , / ) 
NN|TE!6.6)!XMil),RE!l),C0!I)fl«l(Nl) 

6 F0RMAT!5X(FI0.9,10X,Fia.9|l0X,F10.9) 

AU»GAMMA*« I -.b ) 

7 DO 20 I • 1 ,NI 

A • .999*!8./PI )»*.b 
A « 1 . /A 

BKn! I 1 » XM! I )/!RE! I )*AG) 

BKN ! I ) ■ A*BKN ! I ) 

SR » XM ! I 1 *CONVSM 
IF!bR-I.2b) 9b. 9b, 96 
9b COFM » COFMLFISR) 

GO TO 97 

96 COFM ■ C0FMHF!SR) 

97 GM » XM ! I ) 

COCN ■ COCF!GMJ 

DC ■ COFM-COCN ' 

UC ■ CO ! I )-CDCN 
CK ! I ) ■ UC/OC 
WrITE!6,9)0C,UC,CR! I ) 

H FORMAT !/,/, IH ,»DC ■ • , E 20 . 1 0 , b X , ♦ UC ■ • » E 20 . I 0 , b X , • C « ! I) ■*,E20.I0) 
20 CONTINUE 
E » .212 
DO 9 1 J> 1 , 101 
X« .00b* ! J- 1 ) ♦0.3b 
FX*DSUX!CR,RKN,E,N1 ,X) 
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^ x»f X'X 

N 1 WK 1 TE ( 6 , 1 2 I X »NI 
X ■ XF(CH» 

00 M2 J> 1 , lUl 
e ■ .00&* ( J- J ) ♦. I 
Ft>DSOE (CK.BKN|E»Nl,Xi 
^ t«FE-E 

M2 mKI TEI6. 13)E|FE.NI 
X ■ XF ( CM ) 

L > .217 

12 F0RhAT(/,lH ,»X ■» »E20. iO,5X , »FX » • , F 20 • I 0 » 5 X , • N ■ »,I3» 

13 FORMAT!/. IH i*E ■* .£20. 10 iSX . *FE a * . E 20 . 1 0 . 5 X . * N a «,I3| 

S40*n • 0 

DO 7S 1-1 ,N1 
SH-CMaCONVSM 
GMaCM 

IF (SR-1 .2b)as.ub.66 

Bb COLS(I) a CDCF 16M)*| (COFMLF(SR)-COCF!6Fin»EXP(-E/(0KN( 11 )*aXI 1 
GO TO 70 

86 COLS(l) a CDCF (GM)^(C0 FMhF(SW)-C 0CF(6M) )*EXP(-E/(BKN( I ) »aaxi 
70 S^DaSuD«(CDLS( 1 l-CO( n 1«*2. 

7b CONTINUE 

HMSa ( SQD/NI ) .5 

WH 1 TE ( 6 , 7M ) 

7m FORMAT(/,/,/.9X,'COLS’.18X,»CO*,17X.*BXN»,/I 
wHlTE(6,73MCDLS(N),Cn(N),BKN(N).N-l,Nl) 

73 FORMAT! IH ,3E2U. 10) 

>«RITFI6.3)CMtNl .X.E.GSI 

3 FORMAT!/./,/././. IH , 9 X , • C M » , ) B X , * N I • , / , / , 1 H , 7 X , F 6 . M , 1 6 X , I 2 , / . / , / 
I.IM ,9X ,« X •, 19X , »E *, 18X , »6SI •,/,/, IH .3E20.10) 

WH I TE ! 6 , 7 1 )HMS 

71 FORMAT!/,/,/,* R M S - » , E 20 . 1 0 ) 

51 a CM 

52 a CM*C0NV5M 
DCF a CDCF ! 5 1 ) 

DFML a CDFMLF ! 52 ) 

OFMM a CDFMHF!52) 

WKITE!6,76) DCF, OFML, DFHH, 51, 52 
76 FORMAT !/,/,/, IM ,♦ CDCF - • , E 20 . 1 0 . / . / i 1 H ,* C DF ML F » * , E 20 . 1 0 , 

1 /,/,lH ,* CDFMHF-* ,E20. 10,/,/, IH ,* 5 1 ■ * , E 20 , 1 0 , / , / , 1 H , 

I * 52»*,E20.10) 

DO b3 L- 1 ,b 
DO bO Ka 1 , 9 

BKNIX) a ! 0. 0 1 ♦O.O 1 • ! K - 1 ) ) • ! 1 0. •• ! L- 1 ) ) 

BXN!K) a BKN ! K ) a ! 1 0. • • ! -3 ) ) 

5R a CM»C0NV5M 
GM a CM 

IF !5H-1 .25)bb,bb,56 

bb CDLS!K) a CDCF ! GM )♦!! CDFMLF ! SR ) -COCf I n *tXP ! -E/ ! Bkn!K))*«X)) 

GO TO 51 

56 CDLS!K) a CDCF)GMJ4(C0FMHF!5R)-CDCF(6M) )aEXP!-£/( BKN(X))*aX) 

51 WK]TE!6,b2)CM,bR,8KN!K),C0LS!K) 

52 FORMA T ! IHO , 1 6X , *CM* , 1 MX ,* 5H ',/,/, 1 H , I 2 X , F 8 . M , 8 X , F 1 0 , 5 . / , / , / , 1 M , 

1 VX , *BXN* , 1 7X , *COLS* , / ,/ , IH , 3 X , E I 5 , P . 5 X , E 1 5 . 8 ) 

50 continue 

53 CONTINUE 
GO TO 10 
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lOU STOP 

FUNCTION DSUX(CR|HKN,E,NI ,X » 

DIMENSION CR( 100) iBKN( lOU) 

SX « 0,0 
DO 6S I > 1 I N I 
BKNX > BKN ( I ) vtX 

SX ■ SX^(CR( I )-EXP(-E/BXNX ) )*(tXP(-|T/BKNX ) )•( aLOG( 8 < N ( I ) ) ) / B K N X 

6b CONTINUE 

OSUX ■ X*SX 
Pt TUWN 

FUNCTION OSUE<CH#BXN,E ,NI ,X ) 

OlliENSION CR( 100) ,BKN( 100) 

SE > 0.0 
DO 66 I ■ 1 ,N I 
BKNX « BKN ( I ) 99X 

SE ■ SE^(CH( 1 l-EXP(-E/BKNX ) )*(EXP(-E/BKNX ) )/BKNX 

66 continue 

DSUE > E^SE 

RETURN 

END 



PROGRAM AFILIP-LOW CM 


This program computes those same values as the AFILIP-HIGH CM 
except that only Mach numbers lower than unity are employed. 
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WKUN,/TPC AF ILI9.UAHXXXXXXXXX, ORBIT, M,30U 


LOW CM 

OrOR • 1 S MA 1 N • MA 1 N 

XF<S» ■ .399^.016«S 

GSIF(5» ■ .99756 1 7977-. UiaM2696(S*S 

CDFMLF(5»»(2./(PI««.5>)»*<8«/<3.t5)^«.»S/lB.-8.*5*S*5/2l0#)«6Sl 
CJFMHF(S) ■ 2.»( I .♦! ./(S«S)-.25/(S*S*5»S) )*GSI 
COCF(S) ■ 0. 90297-0 . 1 M2M 799*S»S*0. 959^0669* ( • » 

DIMENSION XM(1U0),RE(IU0I,C0(100),BkN(100), 

1 CDLS( lUQ) ,CH( 100) 

PI * 3.19 159265 
GAMMA ■ 1.9 

CONVSM ■ ( G AMM A • . 5 ) • • . 5 

10 RLAOIS, 1 ,END«100)CM,NI 

1 FORMATIFIO.0,15) 

HEA0(5,33) XL.EL.GSl 

33 format ( 3F 1 5.0 ) 

KEA0(5,2> (XM(I),RE(I),C0(1),I«1.N1I 

2 FORMAT ( 3F20.0 ) 

G5 I ■ GS I F ( CM 1 

SELK > 0.0 
SELK2 « 0.0 
SY - j.O 
SYKN b 0.0 
WR I TE ( 6 . 1 1 1 

11 forma T ( I H 1 , 1 3x ,• present experimental sphere oatam 

wR 1 TE I 6 . 5 ) CM 

5 FORMAT! lH0,/,/,25X, *CM - • , F 6 . 9 , / , / . 1 1 X , • X M • , 1 7 X , • R E * , 1 9 X , • C 0 • , / ) 
rRITE(6,6I(XM(I),RC(I).CD(I),I>1,NI) 

6 F0RMAT(5X,F10.9,10X,F10.9,10X,Fi0.9) 

AGbGAMMA*«(-.5) 

7 DO 20 I > 1 ,N I 

A B ,m99b(8./P1 )»».5 
A B 1 ,/A 

BKN ( I ) ■ XM ( I ) / 1 RE ( I ) -AG ) 

BKN ( I I B A«BKN ( I ) 

SR ■ Xrt ( 1 ) »CONVSM 
IF(SR-1.2S) 95,95,96 

95 COFM B COFMLF(SR) 

GU TO 97 

96 CDFM s CDFMHFlSR) 

97 GM B XM ( II 
CDCN B CDCF(GM) 

DC ■ cdfm-cdcn 

OC B CO ( I I -CDCN 

CR ( I I ■ UC/OC 

WR I TE ( 6 , 9 I DC ,0C ,CR ( I ) 

H F ORMa T ( / , / , 1 H ,*0C ■* ,E20. 10,5X , 'UC ■ ♦ , £ 2 0 . 1 0 , 5 X , • C R ( 1 I E20.ini 

20 continue 

t » .212 
DO 91 JB 1 , 1 0 1 
Xb.005* 1 J-1 I ♦0.35 
FXbOSUX(CR,9KN,E,NI ,X) 
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F X«F X-X 

MJ <<NITE(6, 12)X,FX ,M 
X ■ XF ( CM ) 

DO H 2 J ■ I , 0 1 
E«.00^*( J-H*. I 
FE-OSUE ( CM ,8KN ,E ,N I , X ) 

Ft-FE-E 

H2 rtrt I TE ( 6 , I 3 )E tFE , N I 
X ■ XF I CM ) 

E » .212 

12 FORMAT!/, IM ,'X ■ • , E 20 , I 0 , 5 X , • F X » • , E 2 0 . 1 0 , 5 X , • N » »,I3) 

13 FORMAT!/, IM ,»t ■ • , E 20 . 1 0 , 5 X , * F E ■ » , E 2 0 . 1 0 , 5 X , ♦ N > »,I3) 

saD>n,o 
DO 7S I - 1 ,NI 
SR»CM*CONVSM 
CjM»CM 

I f I SR- 1 .25)85,85,86 

85 CDL5!I) « COCF ! GM ) ♦ I ( CDFMLF ! SR ) -COCF ! GM 1 ) •EXP ! -E/ ! 8KN ! I ) 1 ••X ) ) 

GO TO 70 

86 COLSCl) ■ COCF ( GM > ♦ ( CDFMHF ! SR 1 -cncr ( GM ) ) »EXP ( -E/ ! BKN ( 1 ) ) ••X ) 

70 SUO«SQD>!COLS! I )-COl I ) )**2. 

7b continue 

KMS« ! SQD/N I ) •• . 5 
WR 1 TE I 6 , 78 ) 

78 FORMAT!/, /,/,9X, ’COLS’, 18X,'C0»,17X,*BKN»,/) 
mRITE(6,73)!C0LS!N),CD!N),BKN!N),N»i,NI) 

7 3 FORMAT ! 1 H , 3E20. 1 0 ) 

mK 1 TE ! 6 , 3 ) CM ,N I , X ,E ,GS I 

3 FORMAT!/, /,/i/,/»lH ,9X,’CM’,18X,’Nl’». / ,/ ,lH ,7X,F6.8,16X,I2,/,/,/ 

I,1H ,9X ,’ X ’, 1 9X , ’E ♦, 1 8X , ’GSl ’,/,/, IH ,JE20.10) 

I TE ! 6 , 7 1 ) RMS 

71 FORMAT!/,/,/,’ RMS- ’ , E 20 . 1 0 » 

S 1 - CM 

52 ■ CM*CONVSM 
DCF - CDCF !S1 ) 

DFML ■ CDFMLF!52I 
DFMH a C0FMMF1S2) 

R«ITE(6,76) DCF, QFML, DFMH, SI, S2 
76 format !/,/,/» 1 M ,’ CDCF-’ ,F20. 10./»/» IH ,’ COFMLF- ’ , E 20 . 1 0 , 

1 /,/,lH ,’ CDFMhF-’ ,E20. 10,/,/, IH ,’ S1-’ ,E20. 10,/,/, IH , 

I ’ S2-’ ,E20. 10) 

DO 53 L- 1 , 5 
DO 50 K- 1 ,9 

BKN!K) a ! 0.0 1 • 0 1 • ! K- 1 ) ) • ! ! 0. •• I L. 1 ) ) 

RKN!K) a B)(N ! K ) • ! 1 0 . •• ! -3 ) ) 

SR « CM«CONVSM 
GM a CM 

1 F ! SR- 1.25)55,55,56 

5b COLSlK) a COCF ! GM )♦ I ( COFMLF ( SR ) -CDCF ! GM )) *EXP ( -E/ ! BKN!X))««X)) 

GO TO 51 

56 COLSlK) a C0CF!GM)>1C0FMHF!SR)-CDCF(GM))«EXP|-E/! BKN!K))**X) 

51 ..HITE (6,52)CM,SR,BXN{K) ,COLS(K) 

52 F ORMAT ! 1 MO , 1 6X , ’ CM’ , 1 8X , ’ SR ’ , / , / , 1 H , 1 2 X , F 8 . 8 , 8 X , F 1 0 . 5 , / , / , / • 1 H , 

I 9X , ’8KN ’, 1 7X ,’ COLS’ ,/,/, 1 H , 3 X , E 1 5 , 8 , 5 X , E 1 5 . 8 ) 

50 CONTINUE 

53 CONTINUE 
GO TO 10 


125 



100 STOP 

function OSUX(CPtBKN,E*N| ,X I 
DIMENSION CR( lUO) ,BKN( too) 

SX ■ 0.0 
DO 6S I ■ 1 . N I 
BKNX ■ BKN(1)**X 

SX ■ SX ♦ I CR ( I I -E XP ( -E/8KNX ) I • ( E XP ( -E/BKNX n • ( AL06 ( BXN(I)J)/BXNX 

65 CONTINUE 
OSUX > X^SX 

KL turn 

function OSUE (CR i bkn ,e .N I , X ) 
dimension CR( IUO) ,BKN( 100) 

SE ■ 0.0 
DO 66 I>1 »NI 
B^NX > BKN(I)*«X 

SE ■ SE ♦ ( CR ( I ) -E XP ( -E/8KNX ) ) • { EXP I -E/BKNX ) ) /BKNX 

66 continue 
DSUE * E*SE 
RETURN 

END 
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PROGRAM CDCLEV 


This program computes free molecular drag and lift coefficients 
for flat surfaces at various angles of attacks and for a specified 
range of speed ratios. 
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«jKUN,/TP CUCLL»^,UAHXX*XXXAXX,FANTa5TIC-GK,&,1SU 


I'iHLICIT H£AL*b( A-H ,0-Z J 

UldlUSION AN GO (2Q)»CL(2B. 20). 00(2^,20) 

WL AU ( S . I 00 > P J • AL J • X 
S>.2« 1 . J *)2 1 356237 
P I »3 . I M 1 S9265 J 
bwP I ■ I . 7 72H538b09 
b'^2» I . 9 1 •♦2 I 356237 
L*2«71B26ie28959 
AK«X*PI / 180. 

P»0. 32759 1 I 
A I -U. 259829592 
A/«-0. 289996 736 
A3«l .921913/91 
A9.-1 .953152027 
A5» 1.061 H05929 
L0D*9. •3,/PI 
0022 I ■ 1 ,25 
5*= I 

00 22 J« 1 , 1 0 
ANG«X H • ( J- 1 ) 

ANGO ( J I >ANG* 1 80. /P 1 
SAsSiNl ANG) 

CArCOS ( ANG ) 

Z»5*5* .5*SA*SA 
tZ»t •• ( -Z ) 

T»Z/3, 75 
tS»E«* ( -5*5 ) 

«■ 1 .0/ 11 .0*P*S ) 

tHPS«l .O-lAl •R>A2*R*R4A3*1H**3, )^A9«(R»*9. )4A5»1R*»5. ) )*ES 
F»5C/PI«(S*S^1.0-(0.25/lS*5)))*ERFS^(S^10.5/S))»t.S 
C00«2 . •£/ I 5*5*5QP 1 1 
IF (2*3.75) 10. 10. 1 1 

faJ«(Z»«.5)«EZ*(1.0^3,5l66229*T«T ♦3,0899929*(T**9.) 

!♦! .2067992* 1 T*«6, ) ♦ . 2 66 9 7 3 7 • 1 T • • 8 . ) ♦ . 0 36 07 6 8 • ( T • • 1 0 , ) 
2*.0095813*(T**12. ) ) 

H1«(Z**1.5)*EZ*(.5*.8789U599*T*T ♦.51998869*(T**9.) 

2*» 1&089939* 1 T**6. ) ♦ . 0 26 5 8 7 3 3 • 1 T • • 8 , ) ♦ • 00 30 I 5 3 2 • 1 T • • I q • > 
3^. 000329 1 1 * ( T** 1 2. ) ) 

GO TO 12 

eO*. 39899228*. 01 328592/T*. 002253 1 9 / (T*T) 

1- .00157565/(T*T*T)*.00V16281/(T**9,| 

2- . 02057706/1 T**5.)*, 02635537/1 T»*6.) 

3- ,Ul697633/lT**7.)*.00392377/lT**8.) 

B1-, 39899228-. 03988029/T-.00362018/(T*T) 

1*. 00163 80 1/1T*T»T)-.01031555/1T**9.)*.02282967/1T*»5.) 
2-.02«95312/(T**6.)*.01787659/(T**7. )-• 00920059/ 1 T**8. ) 
AL»b, •SA»SA/3. 

«L»-9,«SA*SA*Pl*bA 

C*-PI*.5*5A*lb.*SA*SA/6.) 

A J« AL J 

CLlI,j)»S<i2*SUPI*.5*CA*((BO*Bl)«Aj*(AL*0L»PJ*C*Pj*PJ) 
1*90*1 1, *13. -PJ)*AJ)/(S*S)*01*11.*(1,/3.*PJ/3.)*AJ)/(S*S)) 
C«.(l,J)>(LOD/(l.*LOD))*CL(t,J) 

AJ»- 1 . *9,*Sa*5A/3. 
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i3U*HI-8.*Sa*Sa/3» 

CL»0»3./?.-PMb.*5A«SA/J. 

OU^-.b^Pl/H.-SA^SA/J. 

»C0»AI)4BD*PJ4CC0*PJ*PJ4J0*PJ*PJ«PJ 
1 F Ij- 1 I 20 • 2 I 120 

COf I ,j)»SQ2*SA«( CBQ^Bl )*( l«^Aj*KD)^80»t 1»^AJ*(3.-2.*PJ4.5*PJ«PJ) ) 
l/<2«*S*S)^81*(l.^ALJ*(l./3.*2.*PJ/3,-PJ»PJ/6.))/<2.«S*sn 
?*b«2*(B048l >«(1.4ALJ*(-1»*1*B*PJ*PJ-.5»PJ*PJ»PJ) I/(SA*S*S| 

C0( I , J)«CD( I iJMSQPl 

CO(l,J)»(LOO/(l.^LOO))«CO(I,J)^COO/(l.^LOD) 

GO TO 22 

CJ(I|JJ«SQPM(lt^AJ*J-l.>l.B«PJ*Pj-,B«PJ*PJ*PJ))/S 

CO(|,J)»(LOO/(I.^LOO)»»CO(I,J|»COO/(1.*LOO) 

CONTI NUE 

I TC ( 6 • 200 )PJ t AL J 
MhITE(6.201 MAN60(J)iJ>l||0) 

00 2 I « 1 , 2b 

•~PITE(6t202)I|(CD(l,J),J-1.10) 
mK I TE ( 6 . 203 I PJ » AL J 
WKlTE(6,201MANGO(J),J«ltlO» 

00 3 I * I , 2b 

WRITE(6, 202)1, (CL(I,J),J>I. 10) 

GO TO 50 
FOHMATI 3010.5) 

FOKMATI IMI ,25X,9hC0IS,AN6).8X,3hPJ*,e15*5,HX,MHALJ«,F15.5/) 

F0RrtAT(9X,7HS..ANG«,10£l2.5/) 

forma T ( 3X , 1 2 ,6a , 1 OE 1 2 .5 ) 

FORMATUHl,25X,9MCL(S,ANG).8X,3HPJ»,El5.5,MX,MMALJ«iEl5.5/) 

ENO 
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PROGRAM CFEVAL 


This program computes values of force coefficient for a 
specified speed ratio as a function of gas surface interaction para- 
meters and angles of attack. 
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6HUN,/TP C^£VAL«UAHXXXXXXXXX , F A N T A S T 1 C -GK ,G , ISU 


iHPLlCn wa »L«*»( A-H,0»Z » 

dimension A(S0i2),G(S0t2)|H(S0,2)iQ(S0f2)tP(S0(2) 
dimension E(&Ui2) 

DIMENSION GAMi2G),ALJnS),PJ(lS),CF(20,20,20) 

DIMENSION ASISU,2) •GS(bO,2) .HS(SO,2) iqS(SOi2) .PS(Bn«2) 

DIMENSION FF(bO,2) 

PtAOlS. 100 ISHA 1 
SHA2-SRA 1 *SRA 1 
AN>1 .S780 
AI>t*«M32btMlH6J 
A2* t062606Ql 22U 
A J-.0M7573835H6 
AH«t01 736b06Hbl 
Bl«t2H9983683IU 
B2-.0920UI80037 
B3>t0H069697S26 
bN* .00S2<»N*«9639 
A AU- 1 . 38629H36 1 1 2 
A A I ■0.096663HH2S9 
AA2«0«03b9u0V2383 
A A3-0.037M25637 I 3 
AAN-0.0IN5I 196212 
bbU-U.b 

08 I "0. 1 2M98593S97 
b82*0«0688029bb76 
BB3«0.033283bbJM6 
BBM-O.OONHl 78701 2 
PI-3. mb926b3b9 
0051 1-1,19 

GAMl 1 1-5. 0-( I-l » 

GAMH-5.0*(l-n-P|/ieO.O 
S ING-SIN ( GAMR ) 

C0SA«C05 ( GAMR I 

blHAaSINO 

FMJ >1 ,0-SING-SlNG 

F M2-FM 1 -F M 1 

FM3-FM1-FM2 

FMM-FMl -FH3 

J- 1 

IF (f Ml 121 I 31 |2I 

t < 1 ijl-1 .♦Al-Kh|^A2*FM2^A3-FM3^A‘4*FHM*(ai*FMl^B2-FM2^03-FM3^PM*FMH 
1 ) -ALOi, ( 1 , /F M 1 ) 

FFIl,J»-AA0^AAl-FMl^AA2-FM2*AA3-FM34AAH-FM‘t-lBb0^BBl*FMl^bB2*FM2 
MBB3»FM3 ^BFiH*FM‘ 4)-LOGI I »/FMl ) 

GO TO HI 
a ( 1 ,J)-I . 

F F < I , J » - 1 000. 

A( I,J»-7,4H.-AH*E( 1 ,J)/PI 
Abn,J)-2.>6.*AK-FF(I,J)/PI 

G(l,JJ»AK-M.*E(I,J»/(3.-PI)^AR»16.-cOSA-COSA«FF(I,J|/(9.-Pn 
Gbl I ,J)--AK-FF» I .JI/PI 

H( 1 .J)-M,/3.^Ah-lE( 1 ,J)-(-8R./(9.*Pl )^M.-16.*C0SA*C0SA/(9.*P1 1 ) 
l>PI*SINA*SlNA*.5-8.*FF(I,J)*r0bA-C0sA/(3.*PI)) 
hS( 1 ,J)«M./3.4AR*2.*t( 1 ,J)/(3.?P1 )44R*FF( 1 ,J|«|-1.42./(3.*PI | | 
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1*aK«I-PI»SINA«SINA/H.-PI»3.*(SINA»*«*.)/16.) 

U(liJ>«l./6*«Ah*E(IiJ)*(l6*/(3.*Pl).M.*8«*C0SA*C0SA/(3,*PlM 

l*Ah*(6.*E(I.J»/Pl-«2b*>*I»SlNA»SINA^8.*C0SA*C0SA*FF(I,JJ/<»»«PII) 

ttbUiJ»»l«/6.-AR*E(|.J»*2»/<3.*Pn^AR»FriliJ»*(l.^i6./l3**Pn> 

l-Ah«l.b*FF(I|J)/Pl«>AR«PI*(SlNA*SlNA«(l./B»^3**SlNA«SlNA/32*)) 

PM,J)*-l./6.*AR*EU,JI«(-8./(V,«PlMl.-8.*C05A*C0SA/(9,*Pn) 

I -AR *2 . •£ ( 1 , J ) /R I 

Hb<I.j)*-l./6.^AR*£(I,J)/<3.*Pl|-AR«FF(I|J»*(.25^2t/Pn 

1 ♦ .b*FF ( I , J J *AH/P I 
DO b 1 J« I . 1 1 
DO b 1 K ■ I , 1 1 
AL J ( J ) > • 2* ( J- 1 ) 

Pj(<»«.2*(A-i » 

RJI »PJ < K > 

PJ2-P J 1 -PJI 
P j3»P J I «P J? 

CF(K»Jtn>A(Iil)«AS(I.l)/SRA2«AtD(J)*(G(Iil)«PJl*H(I,l)4PJ2«0ll(|l 

14FJ3*PCI,1M^ALJ<J)«(GS1I,1)4PJ1»HS(I,1)^PJ2»QS(I,1)»PJ3«PS(I,IM 

2/SK A 2 

Cf IK.J.I »>CF(A«J,1 )/( (.♦AR) 
bl CONTINUE 

UO lb J> I , I 1 

MR 1 TE I 6 ,21 1 ) J I AL J( J ) 

211 format I IX.MHALJl • I2,RH) ■ F8.2) 

WRITE! 6,9797 1 
WH 1 TE ( 6 ,9997 ) 

WH I TE 1 6,9899 ) 

DU 1 1 Ib K«1 , 1 1 

mH1TE(6,2111)(CF(K,J,1),1«1,7) 

MHITE(6,2m)(CF(K,J,I),Is8,lH) 
wN|TE( 6.2I1I MCF(K,J,I l,I>lb,19) 

111b WR 1 TE ( 6 , 1 62 I 
IS WK I TE ( 6 • 163 1 
OC TO 99 

162 FORMAT!//) 

163 FORMAT! //// 1 

211* rCKMATCUtl5.9,2XElb.9,2XElS,9,2XElS.9t2XElb,9,2XElb.9,2XElS.91 
97 rORMAT!7x,lHu,16X,lH5,16X,2H|0,lbX,?Mlb,lbX,2M20,lbX,2H2b,lbX,2H3U 
1 1 

9997 F0HMAT(7x,2H3b,lbX,2H90,lbX,2HHb,lbxi2Hb0,lbX, 2Hbb ,lbX,2H60,|bX,2H 
1 6b 1 

9899 F0RMATl7x,2H7U,lSX,2M7b,lbX,2H80,l5x.2M8B,lbX,2H90l 
ICO FORMAT !010.b » 
t NU 
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PROGRAM RUFSPH 


This program was used to compute drag coefficients for a 
number of non-spherically shaped objects. The copy shown here was 
used to compute the drag coefficient of ellipsoids of various eccentric- 
ities and gas surface interaction parameters. 
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yhUN»/TP RUFSPH,UAHXXXXXXXXX,ENViH-0YN-6K,5»l5a 


I^iPLIClT REaL*8U-M,0-Z) 

Oin£NSlON rB<Zl),Ff(2|),THa(2|)tTHF(2l)|TXBC2|),TXF(2l) 
dimension oh 2l I »0R( 21 .21 I 21 I » WI 21 > »'’J<2I » »COF <2I»2l>.<iAM(2|l, 
ITHI (21 ) .MW(2I ) 
dimension ThX(HO) 

OIMCNSION qA( IOI «qA2( IOI 

DIMENSION GR(20t 
DIMENSION F I ( 21 HSS( 21 » 

DIMENSION RH(2n 
X( I I>0*076S26S2| I33SO 
X(2)«0«22778S8St I M US 
X(3U0 i 3737060887|S«42 
X(H)-0.SI086700|9S083 
X IS) ■0*6 3605 3680 7 2652 
XI6)>0*7H633I906H60|5 
X<7I«0*839| 1697182222 
X(e)»0*9l2239««2825l 33 
XI9I-0.96397I9272779I 
XI I0)«0*993|28599185I 
Ml 1 l-o* I527533871307 
M(2)«0* 1991729869726 
WI3)»0* 1920961093183 
M(9)a0*l3|6886389992 
M(5)»0* 1 181995319615 
M(6)«0* 1019301 198172 
Ml 7 ) ■0*0832767915767 
M(B)^0*0626720983391 
M(9)^0*0906019298009 
Ml 10)^0*017619007|32 
P^O* 327591 1 
A I ■O* 259829592 
A2^-0* 289996736 
A3^1 .921913791 
A9^- 1 .953152027 
AS- 1 .061 90S929 
PU3. 1 9 159265358976 
DO 2 1>1 , 10 
J-1 l-I 

gam I I ) -P I • .5* I 1 . -X (J ) ) 

2 MW I I ) (J I 
D031 -1 .10 
J-HIO 

GAM(J)aPl».S*l 1 .♦Xl I ) I 
3 MWIJ ) «w( I I 
DO 90 Mal.s 
ECC ■ O.OOOHM 

GA2(M)-.5*(10***(M-1))*PI/I 360 **3600. ) 

GA I M ) .ECC 
AS- 1 .0 

BS-I 1 .-ECC*ECC)*«.5 
0091 - 1 ,21 
009 J. 1 , 2 I 
9 0R| I , J)-0.0 
DO 29 1-1,10 
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SbA-S|N(GAH( I i I 


CgA>COS(G*H( I ) I 
C2-CGA»CG* 


S2*SGA«SGA 

RH(II«CC2/<BS»BS)»S2/(AS*AS)»»«|-*6) 
0K0TH»((RR(I))'»*3*)*(BS»*«-2»)-AS**J-2.n«CGA«56A 
GN( I |>ATaN(OROTH| 

29 continue 

DO 27 I«l 1 ,20 
GN ( I ) t 0 
NR ( I ) > 1 .0 
27 CONTINUE 
00SI>1 , 10 
SS( I I 
T ■ I . / ( I . ♦p^ss (in 
S*SSI I ) 

EkfS»2.»(S-S«S*S/3**<S«*S.)/I0»*(S»*7,)/M2» 

1/216. I/IPI*». SI 
0 1 ( I )>0.0 
006J« 1 , 20 

Y>ABS(S«C0S(GAM( J) ) ) 

T«l./(I.^P»yl 

ERFY»2.*(T-T»T«Y/3.t(Y«*5.)/l0.-(Y»*7.>/M2. 


{♦(Y**9.)/216.)/IPI**.SI 
CGA>C0S I GAM ( J I ) 

IF (CGA)3l ,33,33 
31 Y*-Y 

ERF Y--ERF Y 
33 PQ»,5^.5*ERFY 

{S«S*Y*(|.4ERFY)45»EXP<-Y«Y>/(PI**,5l 

FIIJ)»IQ»Q»2.*PQ*Q*CGA^PQ»PQ)«»»5 

THI(JI«PI*.5-AC0SI((J»CGA4PQI/FI(J|| 

THX( J)«PI«,5-THI ( J) 

6 continue 
00 70 J> 1 ,20 
ANB-GAM ( J ) -GR ( J ) 

ANF aGAM ( J ) -GR ( J ) 

YbaABS ( S«C0S ( ANB ) | 

YFaAflS(S«COS(ANFI I 
TB»l./( 1 .♦P*YB» 

TFai ,/( 1 .♦PayF ) 

ERFYBa2.«(YB-YB»YB*YB/3.4(YB»«5.)/l0»-|7B*»7.)/92.4(YB*»9 

I (PI*».5» 

tNFYFa2.*(YF-YF*YF*YF/3.*lYF*«&.)/lO»"<YF**7,)/H2.<»(yF*»9 

I IP1**.5) 

C ANBaCOS ( ANB ) 

C ANFaCOS ( ANF ) 

IF (CANS )SI ,S3,53 
Si Y b a .. V B 

EnF YB a-ERF YB 
S3 P(JBa,5^ .S*ERF YB 

«BaS*YB*<l«^ERFYB)^S«EXP(-YB*YBI/(Pl»»»b) 

FB(J|a(uB*0B^2«*PUBaGBaCANB«PQBfP0B)**tS 

THB(J)aPIa,5-AC0SMQBaCANB4PQB)/FB(J)) 

IF (CANF )61 ,63,63 
61 YFa-YF 

ERF YFa-EHF YF 


)/2l6. )/ 
)/2l6. )/ 


135 



63 PWF .5*ERFVF 

QF«S»YF«(l.»ERFYF)4S*EXPj-rF*YF»/<Pl**.5j 
FF(JI»{(JF*QF*2.«P(JF»8 F*CaNF^PQF*PQF)***5 
THF ( JI«PI«*5-AC0S( <UF»CANF*PQF > /ff < J n 
TXBJ J)«PI»,5-THB( J» 

TXF C J ) .p I • .5-THF I J ) 

70 continue 

DO BS L* 1 t20 

0lB«F8a)«C0SlGAM<L)-Pl«*5*THB(L>-GR<L>» 

I*RN|L)*RR(L)*mm(L)«SIN|GAM(U)/COS(gRU)) 

01 I I )«DIB<»OI ( I ) 

8S continue 

01(1 )«0I ( I )*Pl*t5 

007K»1 t21 
PJ(K )•• 1*(K«| ) 

00 8 L- 1 »20 

ThRF«PJ(KI*PM.5*(1#-PJ(X»1«TMF(U 

ORf*FF(L)*COS( TmRF^P l**5-(iAM( L )^GR(l* )*W»ML)*31N 
l(GAM(L»)»RR(L»«RR(LJ/COS(GR(Ln 
0R( 1 ,K )»DR( I ,K )-0RF 

8 continue 

DR( I .K )-Pl«fS*OR( I iK 1 
COF ( I |K )«0R( I |K )/DI ( 1 ) 

7 CONTINUE 
5 CONTINUE 
009|«1 ,10 

MR1TE(6, 101 1 SSI 1) lOI ( 1 ) 

NRITE(6,103) GA(M| 

OOlOKal ,21 

MRITE(6,102) PJ(K),COF(I,K),DR(I,K) 

10 CONTINUE 

9 CONTINUE 

90 continue 

WRITE(6, 105) 

00 28 I»l,20 
AI ■! 

GH( 1 )-GR( I )*180«/PI 

GAM( I )«GAH( I )*180*/P1 

WNITE(6,10M)AI,GAm(I),RR(I),GR(I) 

28 CONTINUE 

105 F0NMAT(//MX,lHl,6X,6HGAH(I),6X,5HRR(I)t6X,5HGR(I)/) 

109 F0RMAT(2X,F10.b,2XiE10«5,2X,E10*5,2xiEl0*5) 

103 format ( 9X , 3HGA- ,E 10*5/ ) 

10 I format (//// 9X»2HS»,ElOt. 2Xi6H01 ( S ) . , E 1 3 . 5/ / 8 x » 2HP J , 6 X , 7HK ( P J , S ) , 

15X,8H0R(PJ,S)/) 

102 F0RMAT(9x,Ft0*5,2X,E10«r . t E 1 0 • 5 ) 

200 F0RMAT(////9X,2HS>«E10«'. ,2X,6H0l(S)«»El0«b,//8X,lHK, lOX ,6HGAM| K ) , 
16X,6HTHI(K),6X,5HFI(K)/) 

201 F URMAT (9X,F 10*5, 2XiE 10*5, 2X,E10#5,2xiEl0«5) 

END 
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PROGRAM CLELAN 


This program was employed to calculate orbit perturbations 
due to lifting satellite shapes. 
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g>rUN,/TP*C CLELANiUAHXXXXXXXXX,0RniT,10i600 


c calculation of SATELLITC orbit change caused by lift and drag 
c perturbations..,., numerical integration using a modified RUNoa- 

c futta technique 

c 

implicit double precision Ia-H,0-Z) 

dimension Yf7),N<7),Q(7l|PHlNT(IOl.PHlN(lO),AAlS)tYY|7t40),0EN(6O 
: ) ,ht( 60) 

double precision mu 
real ESTEP iOIFFEO »R ,0 
external oiffeo 

COMMON Ta ,MU,P ,FF , act ,W,R,aCN|ACK,VMEAN 

C 

C INPUT OF initial VALUES 

C 

c L • number of eccentrictties Times number of perigee heights 

c M ■ number of Different satellite attitqdes 

c N • number of area/mass Ratios 

c NO • upper limit on The number of orbits for a single launch 

READ (5|3) N,M.l.NO 
3 format (MIS) 

READ (S|MMAA(I).I» I.N ) 

M F0RMATIEIS.9) 

Read is, im) iphintik) ,phimik) , k ■ i,mi 
!*♦ format (2FIS.I0) 

REAOIS, 10) VE.RE ,mu 
I 0 format I 3E I S .9 ) 

READ IS, 13) 6A,PJ 
13 format I2F10*B) 

READ (S,|S) I I YYIMMiNN) ,MM • 1,7) ,0EN I NN ) »HT I NN ) , NN ■ l»L) 

IS format <SEIS.9/2EIS.9/2I IS,6) 

PI ■ 3. I M IS926S3S9 
PSTEP-,2 

STEP’X*. 3010299656 
C 

no 1010 I ■ i»N 

AM ■ AA ( I ) 

00 1 020 K ■ I • M 

PI ■ (Pl/IBOt ) • PhINTU) 

P2 ■ IPI/I80. ) • Phi* ir) 

RRITE 16.12) N,MiL,NO 

12 FORMaT(|H|, *FLaT plate in aN elliptic orbit about the earth*, 15X, 

1*416,///) 

RRITE(6|I6) AM.PHIRIF) ,PHINT(F) ,GA,PJ,VE,RE» MU,ESTEP 
16 format 1 1 X , *parameteRs* i /// , * Satellite area/mass ■ » , isx,f7,3,9x 
1 »*M2/KG* .//.♦ attitude ANqLE PHIR ■ • , I 5 X ,F 7 . 3 1 9 X , • DEG * . / / » • ATTIT 
2U0E angle PHINT ■ • , I *♦ X ; F 7 , 3 , 9 X , * DEG • , / / » • GA (GAS SURFACE PARAMET 
3ER) ■ » ,8X ,F8,*4 ,20X , • (GaMMA) X SORT ( I -ALPHA ) ----- REF* I*»//» 
*4* PJ (GAS SURFACE PARAM|TeR) • • , 8 X . F 8 . *4 , // , • EARTm ANGULAR VeLOCI 
sty > • , I I X ,E I *4 . 8 , IX,*RA0/SEC* .//,* EARTH RADIUS ■ •»23X,EI0*6 

6 ,3X . *M* ;// , * Earth gravitational constant ■ • ,7X ,EU *7 ,2X , *H3/SE 

7C2*,//,* eccentric ANOMALY STEP SHE • f » 7X , F I 1 • 7 , 5X » * R AD * * / . I H I ) 
AT«C0S(P2)*C0S(PI ) 

AN»AT*TaN(P1 ) 
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AK«SIN(P2) 


naDaBBBnoBBanaoaaaBnaaBBBanoBaBBaiaaBBBiBiaaaBBaBMBaBaBiDoiBBaD 

main onbital clfments 

VI-ECCENTBIC anomaly 
Y2*SEM|mAjOB AX|S 

Y3-ECCENTB|CITy 

yn.obbital inclination 

TS-ARGUMENT or RERIGEE IRITm RESRECT to ascending node) 
t6»longituoe or ascending node (With resrect to vernal equinox) 
T7*T|ME or RERIGEE RASSaGE 
TIMCbTIME into ORBITIS) 


DO 1000 II - I »L 
DO 20 MM • I ,7 
20 Y(MM) b YYIMM,II> 

RHOO ■ DENI I 1 ) 

HH ■ MT I I 1 ) 

HPI • Y I 2 ) • I I .«Y I 3 ) )-RE 

VPbSQRTI (MU/Y( 2 > )*(I.*Y( 3 |)/(li-Yl 3 m 

OYNARR* RmOO*VR »*2 

OLDYbO.O 

OLDTbO.O 

NORBIT-I 

I TCR -0 

J •« I 

NULL ■ 0 

CALCULATION Or VARIABLES 

2 N NEXTbO 
NOIR • 0 
25 ITERbITER*! 

R ■ Y( 2 )*( I ••Y I 3 )*C 0 SI Yt I n ) 

MbR-RE 

RR«Y( 2 ) • I 1 .-Yl 3 ) ) 

HP ■ RR • RE 

ir (hR*LT. 0 « 0 . 0 R.Y( 3 >«LT* 0 * 0 ) go to RRR 
RHO«RHOO«EXR<“<M"MRI)/MH) 

rr ■ SORT! I 1 .-Ym** 2 I*n •♦YI 3 )*C 0 Sf Y( I I ) )/II •-YC 3 )*C 0 S|Yn » ) ) ) 
R ■ Y( 2 )*( I ••Y( 3 )** 2 ) 

TA ■ ACOSI (COil Yl I ) )-Y» 3 ) )/M ••T| 3 )*C 0 SITI I ) ) ) ) 

TaSIN • ISIN(YII ))»5 0 RT(I.-YI 3 )b* 2 ))/|I*-Y| 3 )*C 0 S(Y| 1 )M ,, 
ir (TASIN) H 0 ( 50,50 
NO Ta«( 2 .*PI )-TA 
50 U«Ta*Y( 5 ) 

ARB||t/rr)*(AT»Y( 3 )*S|N(TA)*AN*||# 4 Y( 3 )*C 0 S(TA))) 

AS«(lf/rr)*(A**II.*Yl 3 )«C 05 (TA))*AN*Y( 3 )*S|NlTA)) 

A«SQRT(AR** 2 *AS** 2 *AX*« 2 ) 

VR« S0RT(HU/Y(2))*(Yl3)fSlN(Y(|)))/(|«-Y(3)*C0SIY(l))) 
VSbS 0 RT(MU*(I**Y( 3 )«C 0 S|Ta))/R)-(VE*R*C 0 S(Y(R))) 
VKbVE*R*COS(U)*cOS(Y(N)I 
VbSQRT IVR** 2 ^VS** 2 ^VXR* 2 ) 
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VMtAN*SO»HMU/r f 
00T«U«*VR4AS*V5*AK*VIC)/(A*V) 

ir (abs<oot);lt,i,oi qo to 51 

DOT* I • 

SI AN6ATK«aC0SI00T ) 

TaS»AN€ATK 
IF fTASI SS«SM,SS 
SM COTAS ■ 0*0 
60 TO S6 

55 COTAS • cotan CTAS) 

56 Taj ■ IFI/2.-TAS) 

CDV«2i*S|N|TAJ)«2.*6A*SlN(TAJ)*C0S((P|/2» ) *P J* C 2 fPJ ) *T a J I 

COR« < VR/V I •COV 

C0S»(VS/V)*C0V 

cok«ivr/v»*cov 

CLZ*-2.*6A*SIN(TAJJ«SlNfPl/2«*PJ*(2.-PJ»*TAJ) 

clR-iar^ivr/vi^sinitajm^clz/sinitasi 

CLS«(AS»»fVS/V)«5lN|TAjn*CL2/S|N(TAS> 

clk«iak«(vk/v»#s!nitaj»i*cli/siwctasi 

CDV ■ CDV* I •RS67T21 I? 

CtZ • COV* ( : I /S1N(P2) I 
AHM ■ AM/2t00 
PL ■ CLZ*AMM 
BO ■ COV*AHM 

ACR*IRhO*V** 2)*|IBL •AR/fVMEAN*SlNCTASM)*U«-T<3)*C0SlYm)l-T|2 
1 )*Y(3)*SINIYn n*IBL ♦ COTaS *80 )/V) 

ACS«(RH0*V**2)*f (BL *A5/ ( VPC AN*S I N ( TaS) )) • ( 1 • -Y O ) *cOS ( Y U I ) ) - 
I CSQRTIi;-YC3)**2»-(VE*C05(Yfm/ VMEANl*l|*-Ym*COSfYlim •*2I* 

I CY(2»/V)*(BL • COTAS *80 I » 

ACR»«RH0*V)*HVF*Y(2»/VmFaNI* SINIYIRM *C0SlUl*m .-Y(3I* cosiyi 
tI}M**2) •(•BL *C0TAS »B0 )♦(BL •V/VHEkN | *AKR ( 1 t»Y < 3 ) *C0S ( Y I) I U ) 
ACN*||,/Fr)*ni.*Yf3»*C0S(TAn*ACR-YmpS|N|TA)*ACS) 
ACT*n./FFI*lTm*8!N|TA)*ACPPn#*Y(3)*C0SlTAn*AC$) 
T|HC*IYm-YI3l*SlN|Y(iy)J/VMEAN ♦YC7) 

RA*Y(2I*( I .♦YI3I ) 

OYNA ■ RM0*V**2 

If HOYN a/OYNAPP* *LT.I *00*») 60 TO SB 

ESTeP • IPSTfPPf OYNA/OY»App|**l*STtftX) »••OI7‘l53^f25l 
60 TO 59 

5S ESTEP** I 396263R0I 
59 0RBIT*N0RB!T 

IF (NDIP.EQ* I •0R.NEXT.E0. 1 ) GO TO 6| 

IF (NORB I T.EO. I • AND. ITER.EQ. I) 60 TO 60 
GO TO 70 
60 OLDRA ■ RA 
OLOT ■ TIME 
TA I -TA 
Y2»Y(2) 

Y3*Y(3) 

YR»y ( M I 
Y5*V(5) 

Y6*Y ( 6 I 
HP I -HP 

oldhp*mp I 

0L0Y3*Y3 
V t aV 
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ACT I «ACT 

write U»tS) (Y(MN) ,HN*I i7) |TA,hR|RhOO 
18 F0RMATI35X, *•••••••*• fWITlAL MAIN ORBITAL CLEMENTS ••••••••• t ,/// 

I I I *. ’ECCENTRIC SEMIMAJOR ECCENTRICITY ORBITAL ARGUMENT 

20f longitude Of*. 8X. ’Time OF’.BX.’TRUE^iSX.’PeRIGEE RERIGCE*»/» 
3* anomaly* | 4X , t AXIS* ,28 X ,• inclination RERIGEE ASCENOINGNOo 
ne perigee Passage anomaly height density • ,/,sx, * Y| # , jox , 
5*Y2*,IIX,»Y3*,I?X,*YH»»I2X,»Y5*,|2X,*Y**.|5X.*Y7*,|2X,*TA*.8 X,»mP* 
6 *7X| *RhOO* t//»FR.4iElR*7« F9«RtFl«»7.2FlR*2.Fl7«7.FlRtS»2clO*N|/./ 
7.52X, * (m-K-S-RaOI ANSI • ;//// ) 

70 IF (Y( I |-(2.*PI ) ) 90.80.80 

80 NORB I T«N0RB I T* I 

81 dclra«oldra*ra 
oelt.time-olot 
0ELTA»2#*P!*I TA»TaI ) 
r>ELA«Y2-Y I 2 ) 
oelact-acti-act 
OCLV.VI-V 
0ELMM»MP*HP I 
0ELC«Ym«Y3 
OELI*Y(*n-YM 
0ELR*Y(5)-Y5 
DCLR»-Y(4|-Y4 
AOECAV ■ OLLRA/TIME 
theta* I 80. *T A/? I 

CCb|80.*Y( I )/PI 

IF ( NOIP.CO. I .ORtNEXTiCQ* II GO TO 420 

OLORA ■ RA 

OLOT ■ TIME 

TA |«TA 

Y2«Y(2) 

Y3-YI3I 

YNbYIS) 

Y5»YI5| 

Y4*Y(4) 

HP I ■HP 

olohp»hp I 

OLOY3-Y3 
V 1 sV 
ACT I -ACT 

IF (NULL.EQ* n 60 TO 82 
Y( I ).Y( I 1-2. •PI 

82 0 YNAP* ( RH0*V**2 I • OOOOI 
0YNAPP»(RH0*V**2) 

OUTPUT 

IF (HP.GT.O.O.ANO.YI 3) .GT.OtO < 60 TO 410 
write (4.700) 

700 format (//,iox,i Satellite has crashed on this orbit or The orbit 
I HAS been reduced to circular*.///) 

410 write (4.800) 

I F (NULL*CO* I ) 60 TO 420 
write (4.801) 

420 WRITE (4.810) NORB I T . THfT A » EE . Y ( 2 ) . Y ( 3 ) . Y ( R ) 
write (4.820) 

WRITE (4.830) Y f S ) . Y ( 4 ) . Y ( 7 ) . T I ME 


141 



r» r» o on 


*W|TE UfSMO) 

«RITC (4|tS0) VRiVS|VK*y,ACNtACT|ACK 
write (6f840> 

write (A|870I R,M,hW.P|WHO.TAS 
write (At880> ADrCAY 

WRiTr U 188 SI DELHP.OCLt lOCLl »DKlW, 0 CLW|| 
WR|TEU» 8 W 0 ) 0 CLTA,DCLA» 0 rLRA iOCLT,DELV,DCLACT,|TER * 

WRITE U.7T5) CLZ»CL>*iCt5,CLR»C0V, CDR » COS . COK . A 1 AR . AS . A« 

IF INULL*E0«n SO TO 1000 
IF (NDIP.EO* I ) 60 TO 2M 

IF (NEXTfEQ* I ) GO TO 2*« 

795 format (37Xt»R»,iSX,*S«(|8X,«K*,/ ,* LIFT C OEFF I C I E NT t , / , 9 X , 90 | 9 • 
no./ .» DRAG coefficient* ,/,9X, 9019. 10./ ,* ATTiTUOE VECTOR*. /. 
|9X .9019, 10.//////) 

800 format l|X.*N0RilT*;i |X,*TA*.17X,*TI*.I7X,*Y2*.l7X.»Y3*,|7X.*y9»,| 
IOX.*TA.YI IN OEgREES*) 

«0I FORMATf|X.*NEXR PERIGEE*) 

RIO format CI5.RX .5019, |2./T 

820 format (I 8X . • YSf . |7X . « Y|* , 1 7X . • Y7* . 14X . «T|ME* ) 

830 format (9X.R0I9. 12./) 

690 FORMAT (eX(*VR*.14X.*VS*.16X.*VK'.UX,*V*.l6X.*ACN<.15X.*ACT*.lSX. 
1 * ACK * ) 

850 format (7018,11,/) 

RAO format (9X,*R*.|9X,*M*,18X,*HP*,I9X,*RA*. 1 8X . *RHO* . 18X . *TaS* ) 

870 FORMAT (AD20*I3.//) 

8R0 format (AX.* the APOGEE DECAY RATE ■ *.1020.T»* METeRS/SECONO • . // ) 
R85 format ( 7X . t oELhP* .OX . tOELE* . 1 RX . •DEL I * . 15X . *OELW» . I 3X . toEUWR* , / , 
I5EI8.I2./ ) 

890 F0RMAT(7X.*0ELTa* .I 3X,*0ELA* ilRX.*OELRA* .1RX.*0ELT*.19X.*0ELV* .OX 
I . *0ELACT* .lOX .•fTCRATlORS* ./.AEI8,12.8X,IS./) 

ITER-0 
J ■ 0 

BnananBDnBBnnaDaBanaBaBBaBBaaaaHBaBBBaaBBBBaanBBBnBBOnaaBBntaanBB 

WO call 6ILL(0IFFE8.Y, , I .ESTeP.W.0.7 ) 

aBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB 

IF (HP,LT,0,0,OR,Y(3) tLEtO.O) GO TO 999 
IF (N0R8IT-N0) 9 I 0 . I 000 . 1 000 

W|0 IF (HP*GT, ( *8*0LDHP ) • ANO« Y( 3 ) ^GT# ( *8*0L0T3) ) SO TC 9)5 
OLDHP-HP 
0LDY3-YI 3) 

NOIP-I 

WRITE (A. 800) 

60 TO 25 

915 IF (0YNA«GT,0YNAP,0R, J.EQ, I ) GO TO 29 
920 Y( 1 )«2.«P;»Y( I ) 

J - I 
NEXT . 1 

WRITE (A.800) 

WRITE (6.802) 

802 format (IX,*NEAR APOGEE*) 

SO TO 25 
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NULL ■ I 

GO To ei 
1000 CONTINUE 
1020 CONTINUE 

1010 continue 
end 

9rOR»l5 OirFEOtOlFFEO 

»EAL function OIFFEO(TiI) 
dimension Y I I ) 

double precision Y ,Ta,mU,R.FF|ACT,U»R|ACNiACK#VMEAN 
COMMON T A, MU, P,FF , act ,U,R .aCN.ACK ,VME AN 
C '4naaHBB0RBnaBBaaaBDaaiiBaaBBBB«BaBaBOOBBaaaBBaBDnBnBaDaBnaBBBBBnao 

C 

C differential equations for orbital parameters, Y 2-Y7, K»|TH RE5PECT 
C To THE independent PARAMETER* TI 

C 

C bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb 

GO to (3iOi32Oi330*3N0iSS0iSA0) I I 
310 DIFFEQ*«2,*CYI2)«*2I/SQRTIMU*P»»*FF*ACT 
GO TO NOO 

320 OIFFEQbISQRTC P/MU)/FF)*(2,*tC0SCTA»*YI3»)*ACT*R*(r.»Y(3)**2)*S|N( 
I Ta ) *ACN/P ) 

GO TO <«00 

330 OIFFEQ»R*COS<U)*ACK/SoRT(MU*P) 

GO TO NOO 

3M0 D|FFE0»IS0RT(P/mUI/Y<3» )• n2.*SlNITAI*ACT»|R/P)*<2,*Yl3UC0SITAl^ 

MY(3)**2)*C0S(TA)I*ACNI/FF*R*Ym»S|N(U)*ACT/(P*TAN(Y(m)n 
IF IY(M)J NOO»3f9,MOO 

350 OlFFEO'IR-SINIUlaACFl/ISORTlP^MUI^SINIYIMn 
IF (YCRn MOO, 399, MOO 

3*0 DIFFE0*fY<2l/HU)*C2.*fP/Yl3)*R*YI3» ) *S I N ( T A ) /FF -3 , bSOR T ( MU/P ) *FF • ( 
IYM»-Ym*SlNIY|m I/VMEAN»*ACT*I (R«PBC0S(TAn/(MU*V(3l*FFn«ACN 
GO TO MOO 
399 D|FFEQ*0,0 
MOO return 
END 

QFOR,|S GILL, GILL 

subroutine gill IDV,Y,2,H,N,Q,N) 

DIMENSION YfNl,RrN),Q(N),AlM),C(R),B(M) 

double precision y 

0ATA(A(|),C(n,BM)*I*l*R)/2**5,2,,2*,2928932B3,l,,2»l*70710*7l, 

II • 1*******6,,5,2,/ 

C 

c This routine is a modified runga-kutta numerical integration 

c technique 

c 
c 
c 
c 
c 
c 

DX-H 
N( I • 

c 
c 
c 
c 


DX. IS The interval size, 

R- IS THE ARRAY USED TO STORE THE 
value of Y*IXI, RMIbFOIXIbI 


FOR THE FIRST INTERVAL THE Q*S ARE SET TO ZERO, FOR 

subsequent intervals The mreviouslt computed q*s 
are used. 
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r% r\ r% 


C - 

DO S J«l tN 
6 0 < J ) aO* 

10 DO 20 J-l 

DO 15 K«2.N 
15 «n()*ovcr ,K-n 

00 20. Ka| tN 

Y<naVllC)*0>t*Afj)*(»lK)-B|J)*0OC)l 
20 0(K)aQ(K)«3«*A(J)»(IMKl»B(J)*Q(K))«C(J)*N(K) 


test if value of INOEFCNOE^T variable 
HAS BECH beached* 

C ......... ................... 

1FIY(1)*DX»GT.2) OXaZ»Y(|) 

IF(DX*LT«ABS(Y( I ) )«2*E«*B) GO TO 25 
!F1Y(1)-Z» |0|25|25 

25 tCTUBN 
END 
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